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Preface 

This  thesis  applies  the  trigonometric  approach  to 
finite  difference  calculus  in  order  to  calculate  column 
buckling  loads.  The  loads  are  found  using  a  virtual  work 
energy  approach.  It  is  hoped  that  my  results  will  encour¬ 
age  the  further  expansion  of  the  trigonometric  method  to 
different  structural  geometries  as  an  alternative  to  finite 
elements . 

I  wish  to  thank  Dr.  Anthony  Palazotto  for  allowing  me 
to  pursue  this  line  of  research  and  for  his  guidance.  In  a 
year  when  he  had  several  thesis  students,  he  always  had  time 
to  help  me.  His  ability  to  provide  information  was  amazing. 

I  also  thank  Captain  Steven  Hannah,  my  predecessor  in 
this  field,  who  provided  guidance  through  his  own  thesis. 
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/  Abstract 

V 

A  trigonometric  approach  to  finite  difference  calculus 
was  applied  to  solve  for  beam  buckling  loads  using  a  virtual 
work  method.  The  trigonometric  equation,  a  truncated  Fourier 
series,  permitted  varying  the  buckling  load  by  adjusting  a 
wavelength  parameter.  Values  for  the  buckling  load  of  a 
variety  of  beams  -  uniform,  homogeneous,  variable  and  dis¬ 
continuous  inertias,  composite  -  were  found  under  a  wide  range 
of  boundary  conditions.  Choosing  the  optimal  wavelength  pro¬ 
duced  the  result  from  literature. An  optimization  scheme  was 
used  which  determines  the  critical  load  by  locating  the  inter¬ 
section  of  two  buckling  load  curves.  The  method  is  accurate 
as  long  as  points  of  interest  -  maximum  inertia,  discontinuit¬ 
ies  -  are  modeled  by  the  nodal  arrangement.  The  trigonometric 
approach  provided  improved  accuracy  over  the  conventional 
approach  for  a  wide  range  of  wavelengths.  For  an  infinite 
value  of  the  wavelength,  the  trigonometric  approach  converges 
to  the  conventional  one.  A  variable  mesh  designed  to  concen¬ 
trate  nodes  about  points  of  interest  was  found  to  be  relative¬ 
ly  ineffective  when  compared  to  a  uniform  mesh.  Composite 
materials  were  modeled  using  an  equivalent  flexural  rigidity. 
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COLUMN  BUCKLING  OF  ISOTROPIC 
AND  COMPOSITE  BEAMS  USING  A 
TRUNCATED  FOURIER  SERIES 


I.  Introduction 


Background 

In  the  past,  structural  problems  were  solved  using  simpli¬ 
fying  assumptions  and  rule  of  thumb.  This  still  holds  true 
today  to  a  great  extent.  But  today's  problems  are  demanding 
more  than  the  "linearized  solution"  often  seen  in  problem  anal¬ 
ysis.  Today's  structures  are  becoming  more  complex,  more  dar¬ 
ing,  more  costly  than  ever  before.  From  new  architecturally 
advanced  domed  stadiums  to  large  space  structures,  these  designs 
push  our  knowledge  and  its  use  to  the  limit.  With  this  advance 
has  come  the  need  for  more  effective  tools  to  provide  the  neces¬ 
sary  accuracy  required  by  these  designs.  And  this  need  is  being 
met  through  the  use  of  modern  computers  and  associated  software. 
This  software  takes  form  as  programs  -  STAGS,  SAP,  NASTRAN  - 
available  to  an  engineer  for  structural  analysis  and  design. 

With  these  new  tools,  solutions  to  the  problems  of  the  real, 
"non-linear"  world  can  be  more  accurately  approached.  The 
only  limit  is  the  precision  of  the  engineer's  model  and  the 
numerical  technique  he  applies.  Because  of  the  increasing 
accuracy  available,  the  engineer  can  concentrate  on  attaining 
the  most  accuracy  for  the  least  time  and  cost  -  symbolized  by 
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choosing  the  best  technique  as  judged  in  these  terms. 

(  Finite  difference  calculus  is  the  technique  applied  in 

this  thesis.  It  models  the  beams  analyzed  within  by  dividing 
them  up  into  a  finite  number  of  degrees  of  freedom  (nodes) 

(5-8) .  A  trigonometric  finite  difference  is  used  as  opposed 
to  a  conventional  finite  difference.  The  conventional  finite 
difference  approximates  a  function's  derivatives  by  combina¬ 
tions  of  a  Taylor  series  expansion  (8-10) .  This  technique  is 
good  when  the  function  in  question  is  polynomial  in  nature. 

The  functions  in  this  thesis  are  directly  dependent  upon  the 
beam's  buckling  shape,  which  is  sinusoidal.  To  model  the 
necessary  derivatives,  a  technique  developed  by  Stein  and 
Housner  (1,2)  and  expanded  upon  by  Hannah  (3,4)  is  incorpor¬ 
ated.  A  truncated  Fourier  series  is  used  to  develop  trigono¬ 
metric  finite  difference  expressions  for  the  needed  derivatives. 
By  using  a  Fourier,  rather  than  a  Taylor  series,  a  closer  ap¬ 
proximation  of  the  derivatives  of  the  sinusoidal  shape  results. 
In  his  thesis,  Hannah,  showed  the  improved  accuracy  of  this 
trigonometric  approach  over  the  conventional  method. 

Beam  buckling,  also  known  as  Euler  column  buckling,  is 
an  important  problem  in  structural  design.  Beam  buckling  is 
referred  to  as  column  buckling  because  the  load  configuration 
on  a  beam  that  causes  axial  buckling  is  most  often  seen  when 
the  beam  is  used  as  a  column,  as  in  a  landing  gear  strut. 

The  prediction  of  the  critical,  or  buckling,  load  on  the  col¬ 
umn  is  of  primary  concern  in  the  beam  design.  For  example, 

1  .  possible  solar  power  satellites,  as  described  by  Nansen  (11), 
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would  be  constructed  of  composite  struts.  With  a  known  load¬ 
ing,  the  problem  in  strut  design  becomes  one  of  finding  the 
optimum  dimensions  to  prevent  column  buckling  while  conserving 
mass  and  space.  The  minimum  force  which  can  cause  buckling 
is  the  critical  buckling  load.  Below  this  load,  a  perfect 
beam  will  remain  flat  in  equilibrium.  At  the  critical  load, 
the  beam  will  deflect  suddenly  to  a  position  of  higher  equilib¬ 
rium.  This  defines  the  buckling  load  as  that  load  which  causes 
the  beam  to  deflect  from  its  lowest  equilibrium  position  (12-13) . 
While  there  are  higher  critical  loads,  they  are  not  of  interest 
in  this  thesis. 

In  Hannah's  work,  the  uniform  beam  was  observed  under  a 
variety  of  boundary  conditions  to  validate  the  use  of  trigono¬ 
metric  finite  difference  expressions.  However,  most  beams  are 
not  so  simple.  Struts  for  use  in  the  solar  satellite  can  be 
both  composite  and  tapered  -  not  constant  in  either  material 
properties  or  inertia.  The  non-uniform  beam  is  of  great  in¬ 
terest.  While  some  solutions  have  been  published,  these  tend 
to  be  expressed  as  bounds  on  the  exact  buckling  load  found  by 
a  numerical  method  or  by  approximations  to  the  governing  dif¬ 
ferential  equations  (14-18).  The  ability  to  achieve  an  accu¬ 
rate  solution  may  be  possible  with  the  trigonometric  technique. 

In  fact,  the  exact  solution  may  be  attainable.  The  solution 
is  only  as  exact  as  the  model  of  beam  buckling  makes  it. 

There  are  many  effects  which  are  not  included  in  the 
beam  buckling  equations.  The  trigonometric  technique  can  pro¬ 
vide  the  solutions  for  beams  with  variable  cross-section, 


material  properties,  node  arrangements  and  combinations  of 
various  boundary  conditions. 

Purpose 

The  purpose  of  this  thesis  is  to  use  the  trigonometric 
finite  difference  approach  to  calculate  the  critical  axial 
load  in  a  beam  under  a  wide  variety  of  conditions.  These  con¬ 
ditions  include  continuous  and  discontinuous  changes  in  inertia, 
isotropic  and  composite  materials,  uniform  and  variable  nodal 
arrangements  and  multiple  boundary  combinations.  Solutions 
are  obtained  using  the  virtual  work  equation  and  are  a  function 
of  a  variable  wavelength  parameter  in  the  truncated  Fourier 
series.  Success  in  applying  this  technique  depends  upon  accu¬ 
racy  of  the  generated  solution  and  the  time  it  takes. 

General  Approach 

In  his  thesis,  Hannah  showed  the  superiority  of  the  virtual 
work  equation  (3,19-21)  versus  the  equilibrium  approach  in  gen¬ 
erating  the  correct  answer  to  the  buckling  load  problem.  Like¬ 
wise,  the  advantage  of  the  truncated  Fourier  series  over  the 
Taylor  series  in  developing  the  finite  difference  derivatives 
for  a  simple  beam  was  demonstrated.  Therefore,  I  started  my 
analysis  of  more  complex  beams  using  the  virtual  work  equation 
and  the  truncated  Fourier  series. 

The  first  and  second  derivatives  of  the  beam's  deflection 
at  an  arbitrary  point  are  developed  from  the  Fourier  series. 

To  integrate  the  virtual  work  equation,  the  beam's  domain  is 
divided  into  a  series  of  nodes,  each  a  constant  distance  apart. 
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The  trapezoidal  rule  is  made  use  of  in  integrating  the  virtual 
work  equation  over  the  length  of  the  beam.  The  solution  de¬ 
pends  upon  the  location  of  the  nodes  and  the  derivatives  pre¬ 
viously  developed.  The  end  result  is  an  eigenvalue  equation, 
with  the  lowest  eigenvalue  related  to  the  critical  buckling 
load.  By  separating  virtual  and  actual  displacements  at  each 
node,  the  eigenvalue  equation  is  transformed  into  an  eigen¬ 
value  related  to  the  buckling  load. 

Imbedded  within  the  algebraic  equations  which  form  the 
eigenvalue  matrix  are  the  effects  of  variable  inertia  and  ma¬ 
terial  properties.  Also,  combinations  of  four  end  restraints 
are  used  in  reducing  the  equations  to  their  final  matrix  form. 
Restraints  considered  are  pinned,  clamped,  guided  and  free. 

A  similar  procedure  is  carried  out  to  examine  the  effect 
of  a  mesh  with  non-constant  intervals  between  node  points. 

This  variable  mesh  requires  the  development  of  finite  differ¬ 
ence  derivatives  which  have  dissimilar  mesh  spacings  on  either 
side  of  the  reference  point.  The  procedure  for  developing 
the  matrix  equation  is  similar  to  that  used  in  the  uniform 
mesh  solution  except  these  new  derivatives  are  inserted  at 
nodes  where  the  mesh  interval  changes. 

Results  pertaining  to  the  eigenvalue  for  a  simple  uniform 
beam  were  obtained  by  Hannah  using  a  method  similar  to  this 
one.  The  solutions  for  his  simple  beam  are  widely  available 
in  literature  (18,19,22-23).  A  broader  range  of  beams  can  be 
analyzed  by  applying  this  thesis.  When  possible,  a  comparison 
with  known  solutions  is  made  (8,14-18,24). 
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II 


Theory 


Assumptions 

The  virtual  work  equation  for  beam  buckling  (19-21)  was 
derived  with  the  following  assumptions: 

1.  Sections  of  the  beam  normal  to  the  longitudinal 
axis  that  are  planar  before  buckling  remain  plane 
during  buckling . 

2.  The  length  of  the  beam  is  much  greater  compared  to 
the  width  or  the  height. 

3.  Displacements  are  small  along  the  beam. 

4.  The  beam  is  flat  before  buckling. 

5.  The  beam  is  composed  of  a  homogeneous,  isotropic 
material. 

A.  For  composite  beams,  interlaminar  effects  are 
ignored  and  the  flexural  rigidity  is  averaged 
through  the  cross-section. 

6.  The  inertia  used  at  a  discontinuity  is  a  weighted 
average  of  inertias  around  the  step  which  preserves 
continuity  and  equilibrium. 

7.  Axial  loading  is  along  the  neutral  axis  with  the 
neutral  axes  of  beam  segments  being  co-linear. 

8.  Only  bending  strain  energy  is  considered. 

The  accuracy  of  the  buckling  load  solution  is  only  as 
good  as  these  assumptions  allow. 
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Virtual  Work  Equation 

The  principle  of  virtual  work  is  used  to  provide  the  basic 
equation  of  this  thesis.  Its  ultimate  solution  is  the  critical 
buckling  load  Pcr>  The  principle  of  virtual  work  states  that 
for  the  variation  of  potential  energy  to  be  zero  (an  equilibrium 
state) ,  the  variation  in  external  work  done  by  the  forces  acting 
on  the  beam  through  a  small  displacement  must  be  equal  and  oppo¬ 
site  to  the  variation  of  internal  work  done  by  internal  forces 
(8) .  Mathematically,  virtual  work  states  that 

6W  =  -<SW.  (2-1) 

e  i 

where  6W  and  SW.  are  the  virtual  external  and  internal  work 
e  x 

due  to  small  displacements.  The  internal  work  can  be  related 
to  a  beam's  strain  energy  by 


6Wi  =  -6U  (2-2) 

The  virtual  work  equation  can  now  be  stated  in  terms  of 
the  external  virtual  work  done  by  the  load  P  and  the  internal 
strain  energy.  For  a  one  dimensional  axially  loaded  beam,  the 
virtual  work  equation  is 


L 

El 

0 


d2w  d26w 
dx2  dx2 


dx 


dw  d6w  , 

1 - T~  dx 

dx  dx 


(2-3) 


The  derivation  of  the  virtual  work  equation  is  shown  in 
Appendix  A 

Figure  1  shows  the  coordinate  system  and  axial  loading 
of  the  beam.  The  beam's  lateral  deflection  is  W,  measured 
from  the  neutral  axis  (generally  the  midplane  for  a  symmetric 
beam).  E  is  the  modulus  of  elasticity,  and  I  is  the  beam's 
inertia. 
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Boundary  Conditions 

Four  boundary  conditions  are  considered  for  beams  in  this 
thesis:  pinned,  clamped,  guided,  free.  These  are  the  most 
commonly  seen  conditions  and  they  can  be  applied  in  any  possi¬ 
ble  combination.  For  the  virtual  work  approach  incorporated 
in  this  thesis,  the  boundary  conditions  used  are 


1. 

Pinned 

W  (A) 

= 

0 

6W(A) 

= 

0 

W"  (A) 

0 

SW"  (A) 

= 

0 

2. 

Clamped 

W  (A) 

0 

SW  (A) 

= 

0 

W'  (A) 

— 

0 

SW’  (A) 

= 

0 

3. 

Guided 

W '  (A) 

= 

0 

5W' (A) 

= 

0 

4. 

Free 

W"  (A) 

= 

0 

SW " (A) 

= 

0 

where  A  is  either  the  boundary  x  =  0  or  x  =  L. 

While  pinned  and  clamped  boundaries  supply  two  conditions 
for  W(A)  or  its  derivatives,  guided  and  free  do  not.  A  third 
order  boundary  condition  on  shear  is  available,  but  because  it 
is  a  third  order  differential  it  is  not  possible  for  use  in  this 
work.  Instead,  an  additional  degree  of  freedom  exists  whenever 
a  guided  or  free  boundary  is  encountered  in  order  to  satisfy  the 
equilibrium  characteristics  at  the  boundary.  This  replaces  the 
shear  condition. 
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III.  Numerical  Technique 
Finite  Difference  Calculus 

The  evaluation  of  the  buckling  load  of  a  beam  can  be  a 
difficult  task,  especially  when  the  variation  of  moment  of 
inertia  or  stiffness  does  not  follow  a  simple  law.  For  these 
reasons,  approximate  methods  of  computation  have  been  devised 
which  provide  bounds  on  the  buckling  load  (6,14-18).  Since 
1910,  a  finite  difference  procedure  for  determining  the  char¬ 
acteristic  values  of  differential  equations  has  been  widely 

investigated.  This  procedure  presents  a  wide  range  of  appli- 

% 

cations  as  well  as  simplicity.  With  finite  differences,  the 
solution  to  beam  buckling  can  be  carried  out  easily,  requiring 
only  the  needed  computing  power. 

Conventional  Approach 

Although  a  truncated  Fourier  series  will  be  used  to  de¬ 
rive  trigonometric  finite  difference  derivatives,  a  study  of 
the  conventional  finite  difference  derivatives  is  useful. 

The  conventional  approach  highlights  many  important  ideas  and 
also  adds  insight  in  applying  the  trigonometric  approach. 

The  concept  of  finite  differences  is  similar  for  both 
the  conventional  and  trigonometric  approaches.  A  function 
under  investigation,  the  beam  deflection  W  in  this  thesis, 
is  modeled  as  a  collection  of  node  points.  These  nodes  can 
be  considered  as  superimposed  upon  the  beam.  Figure  2 
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represents  W(x)  divided  into  a  finite  difference  mesh. 


Fig  2.  Deflection  Shape  Versus  Nodal  Arrangement 
The  nodes  are  numbered,  with  the  deflection  corre¬ 
sponding  to  node  i  as  indicated.  The  mesh  spacing  h  is  the 
distance  between  nodes  and  is  shown  as  a  constant.  Node 
points,  marked  by  the  X's  in  Fig  2,  are  positions  at  which 
the  deflection  is  defined. 

The  conventional  derivatives,  those  derived  using  a  con¬ 
ventional  finite  difference  approach,  are  found  as  combinations 
of  a  Taylor  series.  The  value  of  the  displacement  function  at 
location  X  +  h  can  be  described  by  the  Taylor  series 


W(X  +  h)  -  W(X)  +  hW'(X)  +  jh2W"  (X)  +  |h3W"'(X) 


+  -itf  w^tx)  +  ... 
24 


(3-1) 
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(3-2) 


This  can  be  written  as 

Wi+1  =  wi  +  hwi  +  Jh2wi  +  |b3W?'  +  2Th“w];v  +  ... 

where  i  corresponds  to  the  position  of  node  i  and  +1  is  a 

distance  h  from  i.  Likewise,  -1  is  a  node  a  distance  -h  from 

i.  This  notation  provides 

W(X  -  h)  =  W.  =  W.  -  hW!  +  4h2W?  -  ^h3W*.*‘ 

l-l  i  l  2  i  6  l 

+  2Th4wiV  +  (3_3) 

By  subtracting  Eq  (3-2)  from  Eq  (3-3)  and  rearranging  terms, 


an  expression  for  the  first  derivative  is  obtained. 

wi  "  5h(wi+i  -  -  r  wi  +  • • •  (3-4) 

Typically  all  but  the  first  term  is  ignored,  with  the  second 

term  (h2/6)W7  a  measure  of  the  error  by  truncating.  Adding 

Eq  (3-2)  to  Eq  (3-3)  provides  the  second  derivative  at  node  i 

in  terms  of  the  adjacent  node  points. 

.  W7  =  -j (W  -  2W.  +  W.  .)  -  7~W1V  +  ...  (3-5) 

i  h2  i+l  i  i-l  12  i 

Once  again,  the  first  term  is  the  finite  difference  expression 
with  (h2/12)W^v  representing  the  relative  error  due  to  truncation. 

Because  the  virtual  work  equation  depends  only  on  first  and 
second  derivatives  of  the  deflection  W,  no  further  derivatives 


need  to  be  found.  However,  an  increase  in  accuracy  is  possible. 
The  two  conventional  equations  were  derived  at  a  node  in  terms 
of  deflection  values  at  adjacent  nodes.  However  the  reference 
point  i  can  be  located  between  two  adjacent  nodes.  This  is 
diagrammed  in  Figure  3  and  is  referred  to  as  a  half-station 
approximation . 
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Fig  3.  Beam  Deflection  Using  a  Half-Station  Node  Pattern 


The  development  of  half-station  finite  difference  expres 
sions  follows  the  same  procedure  as  the  full-station  approxi¬ 
mations  by  replacing  ±h  by  ±h/2  and  ±1  by  ±1/2.  The  Taylor 
series  expressions  become 


W(X  +  h/2)  =  Wi+Js  -W  +  |w|  + 


,1V 


+  m)  wr + 


(3-6) 


W(X  -  h/2) 


W.  i. 


W,  - 


k  + 


24V: 


wlv  + 
2J  Wi 


(3-7) 


Subtracting  Eq  (3-6)  from  Eq  (3-7)  provides 


The  error  term  is  represented  by  (h2/24)WV'.  Comparing  with 
the  error  term  from  Eq  (3-4),  (h2/6)WV'f  shows  that  a  first 
derivative  of  W  is  better  represented  by  a  half-station  approx¬ 
imation.  To  produce  an  expression  for  WV,  the  expansions  of 
W(X  +  3/2h)  and  W(X  -  3/2h)  are  needed.  These  prevent  bringing 
into  the  solution.  W^  is  a  reference  point,  not  a  node 
as  was  defined 


“C1  + !)  -  Bi +  &)wi +  Kaf)2wi + 


♦  * 


w 


C1  - 1)  -  -  <rDwi + 1?!)2"!  -  K™)>wl 


(3-9) 


Adding  Eq  (3-9)  to  Eq  (3-10)  gives 


IV 


9h2  27 

W.  ,3  +  W.  3  =  2W.  +  —  WV  +  — W*v  + 
i+|  i-|  l  4  1  64  i 


(3-10) 


(3-11) 


Adding  Eq  (3-6)  and  Eq  (3-7)  provides  an  expression  for  2W^. 
Substituting  into  (3-11)  gives 

1  5h2 


WV  =  — ,(W  ,  -  W  .  -  W.  .  +  W.  ,)  -  —  W1V  +  ...(3-12) 

i  2h2  i+1  i+£  i-f  i-f  24  i 


Comparing  the  error  term  from  this  half-station  derivation 
(5h2/24  W*v)  to  the  error  term  from  the  full-station  derivation 
(h2/12  W^v)  indicated  that  the  full-station  expression  for  WV 
(Eq  (3-5))  is  the  more  accurate. 


The  method  of  producing  derivative  expressions  by  half- 
or  full-station  approximations  can  be  extended  to  trigonometric 


finite  difference  expressions  where  the  comparison  between 
error  terms  is  difficult. 

Trigonometric  Approach 

The  conventional  approach  provides  good  approximations, 
provided  the  function  being  approximated  is  polynomial  in  nature. 
However,  beam  buckling  has  a  sinusoidal  mode  shape.  By  using  a 
Fourier  series,  rather  than  the  conventional  Taylor  series,  one 
should  obtain  better  approximations  due  to  its  trigonometric 
terms.  Stein  and  Housner  (1,2)  recommended  just  this  approach. 
The  recommended  Fourier  series,  a  truncated  one,  was  given  by 

W (X)  =  T  +  T  sin--^  —  Xa-^-  +  t  cos— ^X~X (3-13) 
12  X  3  X 

where  X  is  a  variable  wavelength  parameter  and  X^  is  the  ref¬ 
erence  point,  usually  the  node  i.  The  wavelength  parameter  X 
may  take  on  any  value  from  zero  to  infinity.  Derivation  of 
the  first  and  second  derivatives  follows  closely  from  the  con¬ 
ventional  approaches,  with  much  the  same  terminology. 

Using  the  results  from  the  conventional  approach,  a  half- 
station  approximation  will  be  rsed  for  the  first  derivative  of 
the  deflection  W.  The  derivative  of  Eq  (3-13)  with  respect  to 
X  is 

W'  (X)  =  T2jW(X-X°}  -  T^sin— X^Xq)-  (3-14) 

Evaluating  this  at  the  reference  point  X  =  X  gives 

0 

W' (X  )  =  T  7  (3-15) 

0  2  X 

and 
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(3-16) 


15 


wm i 


i 


a 


Subtracting  Eq  (3-21)  from  Eq  (3-22)  gives  the  expression 
for  Wj 


2X,,  •  irh 
— sin?? 


(Wi+i  -  Wi.i) 


(3-22) 


If  h  =  — sin—,  the  finite  difference  expression  looks  very 

IT  Cm  A 

similar  to  that  derived  by  the  conventional  approach.  This 

h  is  the  equivalent  mesh  size,  as  opposed  to  the  mesh  size  h 

from  the  conventional  approach. 

The  second  derivative  of  the  beam  deflection  W  will  be 

a  full-station  approximation.  Therefore,  the  function  is 

evaluated  at  X  =  X  ±  h.  It  was  shown  earlier  that  the  full- 

o 

station  approach  has  a  smaller  truncation  error.  Evaluating 
W (X)  at  these  points  gives 


W(X  +  h)  =  W.^.  =  W.  +  ^W!sin^ 
o  i+l  l  ir  i  X 


(i)2WJ [l  -  cosl|]  (3-23) 


W(X  -  h)  =  W.  ,  =  W.  -  ^W! sin— 
o  i-1  1  TT  1  X 

+  (?)  l1  -  cost] 

Adding  Eg  (3-23)  to  Eq  (3-24)  provides  W£. 

1 


W"  = 


TThl  '  l+l 


(W  -  2W.  +  W^) 


(3-24) 


(3-25) 


1  t) 

By  using  the  formula  cos20  =  cos20  -  sin20  where  0  =  irh/2X, 
the  quantity 

J(?)2[1-cosIx]  -  (r¥)2sin!?v  (3-26> 


This  revised  expression  is  merely  the  equivalent  mesh  size 
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squared  h2 .  The  second  derivative  can  be  expressed  in  the 


familar  form 


Wi  -  fia  <»i+l  -  2wi  +  "i-l> 


(3-27) 


This  equation  is  similar  to  the  conventional  solution  with  the 
mesh  size  h  replaced  by  the  equivalent  mesh  size  h.  The  sub¬ 
script  i  is  merely  the  reference  point  along  the  beam  at  which 
the  second  derivative  is  required.  For  the  second  derivative, 
the  reference  point  is  a  node.  As  will  be  shown  shortly,  the 
first  derivative  uses  a  reference  point  midway  between  two  nodes. 

Not  only  do  the  trigonometric  finite  differences  greatly 
resemble  the  conventional  derivatives,  they  approach  them  as 
X,  the  variable  wavelength  parameter,  goes  to  infinity.  Stein 
and  Housner  included  X  as  a  parameter  independent  from  the 
process  of  solving  for  the  buckling  load.  It  can  be  chosen 
as  whatever  the  operator  wishes.  The  parameter  does  have  some 
physical  significance,  as  will  be  shown  in  Section  IV.  By 
choosing  different  X's,  multiple  buckling  loads  result,  creat¬ 
ing  hope  that  the  "exact"  solution  may  be  attained  and  not 
merely  bounded.  Choosing  an  optimal  X  will  be  given  further 
attention  later.  The  bound  established  by  the  conventional 
finite  differences  can  easily  be  established.  Using  the  fact 
that  sin0  =  9  for  small  values  of  0,  then  from 


one  obtains 


lim 

h  = 

lim 

2X  .  h 

- Sin^r-r- 

X-*-°° 

X+oo 

TT 

2X 

lim 

h 

lim 

2X 

TTh 

X-*-® 

X-*-» 

it  2X 

(3-28) 

(3-29) 
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Therefore,  for  large  X,  the  equivalent  mesh  h  reduces  to  the 
conventional  mesh  h.  With  this,  the  trigonometric  approxi¬ 
mations  reduce  to  the  conventional  approximations  and  the  low¬ 
er  bound  produced  by  conventional  finite  differences  is  achieved. 
By  knowing  the  conventional  solution,  one  may  obtain  a  range  of 
X  where  the  accuracy  is  better.  Of  added  importance  is  the  use 
of  conventional  procedures  to  determine  the  trigonometric  equa¬ 
tions,  i.e.  use  of  half-station  or  full-station  approximations. 
Because  the  conventional  error  is  known  and  the  two  methods  be¬ 
have  the  same  in  the  limit,  the  same  type  of  error  can  be  ex¬ 
pected  in  the  trigonometric  approach.  This  validates  the  usage 
of  half-  and  full-station  approximations  in  the  first  and  sec¬ 
ond  derivatives. 


Application  to  the  Virtual 
Work  Equation 

The  trigonometric  finite  difference  derivatives  just 
developed  will  be  applied  to  the  virtual  work  equation  derived 
in  Appendix  A: 


dfW  dz5W 
dx2  dx2 


f  d6W 

jf  dx  dx 


dx 


(3-30) 


where  E  is  the  Young's  modulus,  I  is  its  inertia  and  P  is  the 
axial  load  which  will  cause  buckling.  The  equation  is  to  be 
solved  for  the  buckling  load  P  =  P  .  The  virtual  work  equa¬ 
tion  cannot  be  integrated  in  closed  form,  so  numerical  inte¬ 
gration  is  used.  First,  though,  the  beam  is  divided  into  a 
nodal  mesh,  with  the  nodes  numbered  as  shown  in  Figure  4. 
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The  nodes  are  represented  by  X's  along  the  beam  axis. 


Fig  4.  Node  Pattern  for  N  Internal  Nodes 


The  distance  between  nodes  is  the  mesh  size  h  which  equals  the 
beam  length  divided  by  N  +  1.  N  is  the  number  of  internal  nodes 
Nodes  0  and  N  +  1  are  at  the  beam  boundaries.  Nodes  -1  and 
N  +  2  are  fictitious  external  nodes  which  are  required  when 
applying  the  second  derivative  at  a  boundary. 

The  numerical  integration  used  on  the  virtual  work  equa¬ 
tion  is  the  trapezoid  rule.  For  ease  in  presenting  the  solu¬ 
tion,  the  following  definitions  are  made 

fi  *  EiIiW?5Wj  (3-31) 

gJ  -  W!6W!  (3-32) 

i  ii 

where  the  subscript  i  is  the  same  as  used  in  developing  the 
trigonometric  derivatives.  The  virtual  work  equation  can  be 
thought  of  as: 

L  L 

/  fdx  -  P  /  gdx  (3-33) 

0  0 

To  integrate  along  the  beam  length  L,  assume  the  function  f 


is  represented  by  the  curve  in  Figure  5. 


Fig  5.  Numerical  Integration  of  fdk 
Using  Full-Station  Approach 

By  using  the  trapezoidal  rule  on  Figure  5,  the  integral  of 


fdX  can  be  represented  as  a  summation, over  the  beam's  nodes: 


Likewise,  numerical  integration  can  be  applied  to  the  integral 
of  gdx.  This  integration  depends  upon  the  values  of  g  at 
points  half-way  between  the  node  points  because  g  is  a  func¬ 
tion  of  half-station  approximations.  Figure  6  demonstrates 
this  integration. 
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g 


Fig  6.  Numerical  Integration  of  gdx 
Using  Half-Station  Approach 
The  integral  can  be  represented  by 


N 

h  1  gi+A 
i*0  1  2 

Replacing  the  virtual  work  integrals  by  these  new  expressions 
gives  the  equation 


h(K  +  t  fi  +  H*l)  '  Ph  I  9i+* 

i*l  i-0 


(3-34) 


The  mesh  size  h  is  factored  out  of  both  sides  of  Eq  (3-34) . 
Using  Eq  (3-31)  and  Eq  (3-32)  and  substituting  the  trigono¬ 
metric  derivatives,  a  definition  for  f  ,  g  is  produced 

Jv  Jv 
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E  1 

f.  =  -4-^  (W,  -2W  +  W,  , )  ( 6W  -  2SW  +  «$W,  J  (3-35) 


k+1 


k-1 


k+1 


k-1' 


9k  -  Si  °w  -  V*»  “V  -  SV±> 


(3-36) 


*  JV-r  2  rw  2  rv  I  2  rw  2 

for  k  =  0  to  N  +  1.  The  expressions  in  Eq  (3-34)  may  now  be 
defined  in  terms  of  the  beam  displacement  W,  stiffness  E,  and 
inertia  I  at  each  node  in  the  beam.  The  solution  of  this  equa¬ 
tion  will  be  the  critical  buckling  load. 

The  summations  in  Eq  (3-34)  can  be  expressed  using  Eq 
(3-35)  and  Eq  (3-36).  By  regrouping  the  coefficients  of 
specific  virtual  displacements  6VT,  the  equation 

a  (P,E,I,W  ,  W  ,W  )<SW  +  a  (P,E, I,W  ,W  ,W  ,w  )6W  +  . 

-1  -101-1  0  -10120 

...  +  a  (P,E,I,W  ,  W  ,  W  .W  )6W  , 

N+l  N-l  N  N+l  N+2  N+l 


taN+2(P'  WN'Vl'V2>%+2  -  0 


(3-37) 


is  produced.  The  a^  are  functions  of  the  axial  load  P  and  the 
actual  displacements  W  in  the  neighborhood  of  node  i.  Also, 
the  a^  are  functions  of  E  and  I  in  the  vicinity  of  i.  Eq 
(3-37)  can  be  written  in  the  form 

a.6W.  =  0  (3-38) 

i  l 

1 


The  trivial  solution  =  0  is  ignored.  The  internal  virtual 
displacements  are  arbitrary  and  are  taken  as  non-zero.  At  the 
beam  ends,  the  virtual  work  idea  makes  actual  and  virtual  dis¬ 
placements  compatible.  For  example,  a  beam  pinned  at  the  left 
boundary  has  the  real  boundary  conditions  W(0)  =  0  and 
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W"(0)  =0.  Because  virtual  boundary  displacements  must  be 

compatible,  6W(0)  =  0  and  <$W"(Q)  =  0.  Applying  the  finite 

difference  derivatives  gives 

<5W  =  0  (3-39) 

0 

<SW  =  -<SW  (3-40) 

- 1 

By  applying  the  boundary  conditions,  external  virtual  displace¬ 
ments  can  be  related  to  internal  and  boundary  virtual  displace¬ 
ments.  Virtual  displacements  at  the  boundary  are  zero  if  pin¬ 
ned  or  clamped  and  undefined  if  free  or  guided.  In  this  latter 

case,  the  virtual  boundary  displacement  <$W  or  SW  ,  depending 

o  N+l 

on  which  boundary,  produces  an  additional  unknown  and  an  addi¬ 
tional  equation.  Applying  the  boundary  equations  reduces  Eq 
(3-38)  to  an  equation  with  N  coefficients  (N+l  with  a  free  or 
guided  boundary,  N  +  2  with  both  boundaries  free  or  guided) 

N 

£  =  0  (3-41) 

i=l 

The  b^  are  functions  of  W,  P,  E,  I  as  were  the  a^  in  Eq  (3-38) . 
However,  the  actual  boundary  conditions  are  applied  to  the 
actual  displacements  in  the  a^  in  the  same  manner  as  the  vir¬ 
tual  displacements.  For  a  pinned  left  boundary,  the  conditions 
are 

W  =  0  (3-42) 

o 

W  =  -W  (3-43) 

”i  l 

Applying  the  boundary  equations,  the  expressions  a^  are  re¬ 
duced  to  b^.  Once  again,  a  free  or  guided  boundary  displace¬ 
ment  is  left  as  an  unknown. 

In  Eq  (3-41),  all  the  virtual  displacement  may  now  be 


treated  as  arbitrary.  Because  of  this,  the  coefficients 
must  be  zero.  This  reduces  the  virtual  work  equation  to  a 
set  of  N  equations  (N  +  q  where  q  is  the  number  of  free/ 
guided  boundaries)  of  the  form 

b±  =  0  i  =  1,2,.. .,N  (3-44) 

The  b^  are  functions  of  P,  E,  I  and  the  actual  node  displace¬ 
ments  along- the  beam.  By  factoring  out  the  actual  displace¬ 
ments,  Eq  (3-44;  can  be  represented  in  its  matrix  form. 


The  elements  of  C  depend  only  upon  E,  I  and  P.  Matrix  C  is 
of  size  N  x  N  and  is  populated  along  and  near  to  the  diagonal. 
Once  the  values  of  E  and  I  are  known  for  each  node  in  the  beam, 
the  eigenvalues  of  the  matrix  C  can  be  solved  for.  Because 
Eq  (3-45)  is  a  set  of  homogeneous  equations,  setting  the  de¬ 
terminant  of  C  equal  to  zero  will  produce  N  eigenvalues.  Al¬ 
though  the  solution  for  buckling  load  yields  an  infinite  number 
of  values,  only  N  are  obtained  from  the  determinant.  This  is 
due  to  the  approximation  of  the  beam  continuum  by  N  internal 
node  points.  By  ordering  the  eigenvalues,  the  critical  buck¬ 
ling  load  is  found.  It  is  P  ,  the  lowest  non-zero  eigenvalue. 
A  zero  eigenvalue  corresponds  to  rigid  body  motion,  a  non¬ 
buckling  phenomenon. 
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Now  that  the  C  matrix  is  determined,  it  merely  awaits 
some  type  of  reduction  to  yield  the  buckling  load.  One 
simple  technique  involves  checking  the  sign  of  the  determinant 
of  C.  A  guess  is  made  at  .the  value  of  the  axial  load  P.  It 
is  inserted  into  the  determinant  and  the  sign  of  the  determi¬ 
nant  of  C  is  checked.  A  larger  P  is  assumed  and  the  process 
is  repeated.  When  the  sign  of  the  determinant  of  C  changes, 
two  bounds,  one  high  and  one  low,  are  established  for  P.  This 
process  repeats  until  the  determinant  of  C  equals  zero.  The 
value  of  P  for  this  condition  is  an  eigenvalue  of  Eq  (3-45) . 

If  it  is  the  lowest  eigenvalue,  then  P  is  the  critical  buck¬ 
ling  load  P  .  However,  as  Franklin  (26)  shows,  there  are 

C  3T 

better  ways  to  solve  for  the  eigenvalues  of  the  C  matrix. 

This  simple  algorithm,  though,  helps  to  illustrate  how  it 
could  be  done.  More  sophisticated  techniques  for  solving  the 
eigenvalue  problem  include  the  method  of  elimination  by  tri- 
angularization,  the  method  of  Householder  and  Bauer,  and  the 
LR  and  QR  methods  (25,26)  . 

The  method  used  in  determining  both  the  critical  buckling 
load  (eigenvalue)  and  the  corresponding  mode  shape  (eigenvector) 
is  from  a  prepared  routine  available  in  the  International 
Mathematical  and  Statistical  Libraries  (IMSL) .  To  use  this 
routine,  Eq  (3-41)  must  be  rewritten.  The  C  matrix  in  Eq 
(3-41)  contained  E,  I  and  the  critical  force  P.  This  force  P 
was  evident  in  the  right  hand  side  of  the  virtual  work  equation 
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(3-30) .  Because  it  appears  on  only  one  side  of  the  equality, 
the  P  dependence  can  be  factored  out  of  the  C  matrix  into  the 
form 


=  P  I  B 


(3-46) 


The  A  matrix  corresponds  to  /  El  W"  6W"dx  and  is  dependent 

upon  E,I.  The  B  matrix  corresponds  iO  the  right  side  of  the 

L 

virtual  work  equation  P  ^  W'SW'dx  and  depends  only  upon  P 
which  is  factored  out  of  the  equation.  It  is  in  the  form  of 
Eq  (3-46)  that  the  prepared  IMSL  routine  solves  for  the  eigen¬ 
value.  IMSL  does  so  by  first  reducing  A  to  upper  Hassenberg 


form  and  B  to  an  upper  triangular  form.  A  is  then  further 


transformed  to  a  quasi-upper  triangular  form  which  is  an  upper 
Hassenberg  form  with  no  two  consecutive  subdiagonal  elements 
being  zero  (31) .  In  this  form,  the  eigenvalues  are  returned. 

In  addition  to  the  eigenvalue  Pcr,  the  corresponding 
eigenvector  W  is  returned  from  the  IMSL  routine.  The  eigen¬ 
vector  is  found  by  using  an  extension  of  LR  and  QR  triangular- 
izations.  The  returned  eigenvector  describes  the  displacement 
at  each  internal  node  of  the  beam  as  well  as  those  at  guided 
or  free  boundaries.  The  displacements  are  normalized  with  the 
maximum  displacement  defined  as  1.0.  The  displacement  curve 
is  available  as  a  by-product  of  the  eigenvalue  analysis.  By 
using  a  prepared  curve  fitting  program,  the  displacement 
curve  is  drawn  through  the  known  node  displacements.  This 
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provides  a  method  for  looking  at  the  effect  a  change  in  the 
beam  has  on  its  mode  shape. 


In  the  guided-free,  guided-guided  and  free-free  beams, 
rigid  body  motion  can  show  up  in  both  eigenvalue  and  eigen¬ 
vector  analysis.  In  the  eigenvalue  problem,  a  critical  buck¬ 
ling  load  of  zero  is  returned  as  a  result  of  the  right  body 
translation.  ‘The  computer  is  merely  programmed  to  avoid  this 
zero  and  return  the  next  lowest  eigenvalue.  In  eigenvector 
plots,  rigid  body  motion  affects  the  displacement  curve  by 
introducing  a  scaling  factor.  To  avoid  this  in  the  guided- 
free  beam,  all  displacements  are  re-scaled  by  subtracting  the 
minimum  nodal  displacement  (usually  at  the  guided  boundary) 
from  a  nodal  displacement  value  and  dividing  this  difference 
by  (1  -  Wm^n  ) .  This  rescales  the  nodal  displacements  and 
produces  the  expected  curve.  A  similar  technique  is  used  for 
the  free-free  and  guided-guided  beam.  An  additional  change 
is  needed  because  the  center  of  the  beam,  at  x  =  L/2,  does 
not  have  zero  transverse  displacement.  This  arbitrary  trans¬ 
lation  is  removed  from  the  nodal  displacements,  producing 
the  expected  curve. 

Summary 

The  eigenvalue  matrix  problem  is  now  set  up.  It  was  de¬ 
rived  from  the  use  of  a  truncated  Fourier  series,  featuring 
the  variable  wavelength  parameter  X,  in  developing  finite 
difference  expressions.  These  expressions  were  applied  to  a 
beam  at  several  nodes,  permitting  the  numerical  integration 
of  the  virtual  work  equation.  The  equations  (3-41)  and  (3-42) 


are  the  end  result  of  this  integration.  By  solving  for  the 
determinant  of  the  C  matrix  in  Eq  (3-45)  ,  the  eigenvalue  P 

Gi 

and  eigenvector  W  are  obtained. 

The  forth  section  of  the  thesis  will  apply  this  technique 
to  simple,  uniform  beams  to  verify  the  method  before  looking 
at  more  difficult  problems.  Besides  verifying  the  approach, 
an  attempt  to  determine  a  best  wavelength  parameter  X  will 
begin. 
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IV.  Uniform  Beam  Results 

The  first  beams  that  were  examined  were  the  simple,  uni¬ 
form  beams.  These  are  beams  which  have  a  constant  stiffness 
E  and  a  constant  moment  of  inertia  I.  Equation  (3-46)  was 
programmed  for  a  computer.  It  allows  for  N  internal  nodes 
with  constant  or  varying  E  and  I.  To  test  the  program,  it 
was  run  for  the  uniform  beam  as  described  above.  Because  E 
and  I  are  constant,  only  the  boundary  conditions  are  different 
from  beam  to  beam.  The  10  boundary  conditions  considered  are: 

1.  Pinned-pinned 

2.  Pinned-clamped 

3.  Pinned-guided 

4.  Pinned-free 

5.  Clamped-clamped 

6.  Clamped-guided 

7.  Clamped-free 

8.  Guided-guided 

9.  Guided-free 

10.  Free-free 

For  these  10  cases,  the  exact  critical  buckling  load  and  mode 
shape  are  available  through  solving  a  governing  differential 
equation  which  was  not  applied  in  this  thesis.  These  results 
are  well  catalogued  in  literature  (18-20,23,24). 

Accuracy  of  the  computer  program  has  been  checked  by 
comparing  the  calculated  buckling  load  and  mode  shape  with 
the  exact  solutions.  As  discussed  in  Section  III,  the  mode 
shape  is  actually  a  curve  fitted  through  the  deflection  at 
each  node  of  the  beam.  This  mode  shape  is  produced  indepen¬ 
dent  of  the  wavelength  parameter  X.  To  show  this,  a  discussion 
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of  the  accuracy  of  the  critical  buckling  load  must  first  be 
made. 

Eq  (3-46)  consists  of  two  matrices  A  and  B.  They  were 

formed  by  numerically  integrating  the  virtual  work  equation. 

Imbedded  within  both  matrices  is  a  dependence  Upon  W  ,  W  , 

12 

. ..,  and  the  corresponding  virtual  displacements.  These 
were  factored  out,  leaving  Eq  (3-46) .  The  A  matrix,  consist¬ 
ing  of  f^  (Eq  (3-35))  applications,  now  depends  upon  E^I^  and 
l/h4.  Likewise,  the  B  matrix  depends  upon  multiples  of  1/h2, 
obtained  by  applying  Eq  (3-36) .  The  equivalent  mesh  size 

<■«  2  X  TTh 

h  =  — sin-^— .  For  the  uniform  beam,  solutions  are  obtained 
with  a  "uniform"  mesh.  By  choosing  the  desired  number  of  nodes 
N  within  the  beam,  the  beam  is  divided  into  a  mesh  with  the 
distance  between  each  node  equal  to  h,  where  h  =  (beam  length)/ 
(N  +  1) .  Therefore,  for  a  uniform  mesh,  h  is  constant  through¬ 
out  the  beam.  This  was  used  earlier  when  h  was  divided  out  of 
Eq  (3-34).  As  will  be  shown  in  Section  VI,  this  simplifica¬ 
tion  is  not  permitted  when  the  mesh  varies  in  size  throughout 
the  beam.  Using  the  fact  that  h  is  a  constant,  fixed  in  size 
when  N  is  chosen,  h  is  dependent  only  upon  the  mesh  size  and 
is  constant  for  a  fixed  wavelength.  The  equivalent  mesh  is 
only  a  function  of  h  and  X,  which  are  determined  by  the  opera¬ 
tor.  The  wavelength  parameter  is  independent  of  E,  I. 

With  the  simplification  of  the  uniform  mesh  beam,  l/h4 
and  l/h2  will  be  factored  out  of  the  A  and  B  matrices  of  Eq 
(3-46) .  This  produces  two  matrices  A  and  B. 
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results  in  an  equation  similar  in  form  to  Eq  (3-46) . 


fc 


This  equation  highlights  one  of  the  main  advantages  of  Stein 
and  Housner's  truncated  Fourier  series.  The  eigenvalue  equa¬ 
tion  does  not  return  a  single  critical  buckling  load  for  N 
nodes  as  does  the  conventional  finite  difference  equation. 
Instead,  the  eigenvalue  equation  returns  a  value  for  kcr  which 
permits  Pcr  to  become  a  function  of  X.  By  choosing  a  X,  a 
critical  load  is  returned.  This  is  the  basic  advantage  of  the 
Fourier  series  coupled  with  X.  By  choosing  a  good  X  it  may  be 

possible  to  obtain  a  P  better  (more  accurate)  than  the  con- 

cr 

ventional  finite  difference  method  generates.  And  if  X  is  just 
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right,  the  "exact"  solution  may  be  obtained.  The  question  of 
what  is  a  good  wavelength  parameter  will  be  discussed  further 
in  this  Section. 

In  addition  to  determining  the  buckling  load,  it  is  some¬ 
times  desirable  to  determine  the  corresponding  mode  shape.  By 
using  the  IMSL  routine  described  in  Section  III,  the  eigen¬ 
vector  can  be  determined  for  either  Eq  (3-46)  or  Eq  (4-1) . 

For  the  uniform  mesh  beam,  Eq  (4-1)  is  used.  The  eigenvector 
corresponding  to  each  of  the  N  eigenvalues  is  stored  and  is 
available  upon  request.  Because  this  thesis  is  mainly  con¬ 
cerned  with  eigenvalue  determination,  the  methods  of  eigen¬ 
vector  determination  will  not  be  discussed.  For  my  purpose, 
if  an  eigenvector  is  needed,  it  will  be  that  one  returned  by 
the  IMSL  package.  The  eigenvector  is  of  interest  though, 
because  it  is  independent  of  the  wavelength  parameter.  As 
was  just  shown,  the  A  and  B  matrices  have  no  direct  X  depen¬ 
dence.  Instead,  they  depend  on  the  buckling  parameter  kcr 
which  was  a  constant  in  the  eigenvalue  solution.  Therefore, 

•  the  eigenvector  delivered  depends  upon  A,  B  and  kcr.  The 
eigenvector  is  independent  of  X.  In  the  computer  program, 
the  buckling  mode  shape  is  determined  before  the  operator 
decides  on  what  X  or  range  thereof  he  is  interested  in.  Be¬ 
cause  of  this,  even  though  a  poor  X  may  be  chosen,  resulting 
in  a  poor  Pcr#  the  mode  shape  is  independent  of  the  X. 

Uniform  Beam  Solutions 

As  stated  in  the  introduction  to  this  section,  10  cases 
of  a  uniform  beam  were  studied  to  verify  the  designed  computer 
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program.  These  10  cases  differed  only  in  the  boundary  condi¬ 
tions,  the  wavelength  parameter  X  was  usually  varied  between 
0.2  and  2.0.  This  produces  a  range  of  P  ,  according  to  Eq 
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(4-3) .  Also,  the  conventional  finite  difference  solution  was 
calculated  to  provide  a  comparison  with  the  trigonometric  sol¬ 
ution.  The  conventional  solution  can  be  found  by  letting  X-*®, 
or  in  the  uniform  beam  case,  replacing  h  by  h.  Because  h  was 
factored  from  the  matrix  equation  (4-1) ,  this  substitution  is 
easily  accomplished. 

Figures  9  through  12  at  the  end  of  this  Section  compare 
the  buckling  load  curve,  produced  by  varying  X,  with  the  con¬ 
ventional  solution.  Because  the  correct  solution  is  known, 
the  percent  error  of  the  buckling  load  is  plotted  versus  the 
wavelength  parameter.  Because  the  conventional  finite  differ¬ 
ence  solution  is  independent  of  X,  it  is  represented  as  a 
dashed  line  of  constant  error.  Also,  a  zero  error  line  is 
included  which  allows  for  ease  in  comparison  and  in  determining 
the  best  wavelength  parameter.  Although  10  cases  were  examined, 
several  combinations  of  boundary  conditions  produced  the  same 
curve.  As  an  example,  pinned-pinned  and  free-free  boundaries 
produce  the  same  buckling  load  curve.  Therefore,  only  the 
four  independent  curves  are  included.  The  boundary  combina¬ 
tions  they  represent  are  included  in  the  figure  title. 

In  addition  to  buckling  load  curves,  the  calculated 
eigenvector  is  plotted  for  five  boundary  combinations.  Once 
again,  one  eigenvector  can  represent  several  different  boundary 
combinations.  As  a  check  on  the  accuracy  of  the  eigenvector. 
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the  exact  buckling  shapes  are  plotted  on  the  same  figure.  Re¬ 
sults  are  shown  in  Figures  13  through  18. 

For  ease  in  calculation,  several  parameters  were  normal¬ 
ized.  For  simplicity,  E  and  I,  variables  within  the  A  matrix, 
were  defined  as  1.0.  Because  the  beam  is  uniform,  E  and  I  do 
not  change.  They  are  represented  as  the  ratios  of  E/Emax  and 
I/Imax  which  accounts  for  the  normalization. 

The  beam  length  L  is  also  normalized  as  L/Lt)eain  “  1.0. 

In  Section  V,  discontinuities  are  introduced  which  make  this 
practice  of  normalization  quite  useful.  Also,  the  user  is  not 
forced  to  choose  specific  numerical  values  for  E,  I  and  L. 

The  critical  buckling  load  is  represented  by  Pcr  =  m  EI/L2 
with  the  value  of  m  returned  by  Eq  (4-3) .  Once  the  solution 
is  in  this  form,  the  user  can  choose  any  values  of  E,  I  and  L 
desired,  without  having  to  solve  Eq  (4-2)  for  each  specific 
case. 

Using  these  techniques,  the  10  uniform  cases  were  solved 
for  Pcr  and  the  buckling  mode  shape.  A  few  of  the  results 
are  discussed  here. 

Pinned-Pinned 

The  most  common  column  buckling  problem  is  the  uniform 
pinned-pinned  column.  This  problem  is  usually  the  first  buck¬ 
ling  problem  a  student  sees.  For  analysis,  two  cases  were  run. 
Nine  internal  node  points  were  used  in  one  case.  Figure  9,  and 
19  in  another.  Choosing  these  numbers  of  internal  nodes  re¬ 
sults  in  a  mesh  size  of  0.1  and  0.05,  for  a  beam  of  length 
equal  to  one.  It  can  be  seen  that  19  nodes  provides  a  mesh 
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twice  as  fine  as  for  N  -  9  nodes.  The  nine  node  beam  is 


shown  in  Figure  7  with  its  buckled  mode  shape  superimposed. 


formula 

PT 

P  *  (mr)2ry  n  =  1,2,3,...  (4-4) 

cr  Ju 

The  lowest  Pcr  is  given  for  n  =  1  which  produces  Pcr  =  tt2  EI/L2 

*  9.8696044  EI/L2.  The  nine  node  solution  results  in  values  of 
Pcr  more  accurate  than  the  conventional  finite  difference  for 
the  range  of  .71  <  X  <  •,  approximately.  Also,  the  exact  buck¬ 
ling  load  is  calculated  when  \  *  1.0.  At  this  value  of  the 

* 

wavelength  parameter,  Eq  (4-3)  returns  m  =  9.8696044  which 
matches  the  known  solution  for  n  =  1. 

In  addition  to  matching  the  buckling  load,  the  eigenvector 
that  is  produced  exactly  matches  the  known  mode  shape  W(X/L) 

*  sin  ir(X/L).  Actually,  the  nodal  displacements  lie  on  this 
curve,  with  the  computer  producing  a  curve  through  them. 

For  19  internal  nodes  the  same  results  are  obtained,  but 
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with  greater  accuracy.  Once  more,  at  X  =  1.0  the  exact  solution 
is  returned.  For  the  range  . 71  <  X  <  ®,  Pcr  is  more  accurate 
than  the  conventional  solution.  For  this  same  range,  the  19 
node  buckling  curve  is  more  accurate  than  the  nine  node  curve. 
Finally,  the  eigenvector  matches  the  known  mode  shape  again. 

In  both  cases,  X  =  1.0  resulted  in  the  exact  solution. 

This  value  of  X  can  be  related  to  the  known  buckled  shape.  If 
the  beam  problem  were  solved  explicitly  with  L  non-dimensional- 
ized,  the  exact  solution  would  occur  at  X  =  1.0L.  The  fact  that 
X  is  a  wavelength  parameter  gives  a  clue  to  its  nature.  If  the 
sine  curve  of  the  buckling  mode  shape  were  extended  beyond  X  =  L, 
it  would  have  a  wavelength  of  2L.  Its  half -buckled  wavelength 
would  be  L,  which  corresponds  to  the  optimum  value  of  X.  The 
idea  of  XQpt  corresponding  to  the  half-buckled  wavelength  will 
be  checked  further. 

Guided-Guided 

A  second  type  of  beam  that  was  checked  is  one  in  which 
both  boundaries  are  guided,  that  is,  the  boundary  is  free  to 
translate  vertically,  but  must  retain  zero  slope  (dw/dx)  at 
each  end.  Although  the  beam  is  not  common,  it  does  provide  a 
check  on  solutions  produced  for  more  exotic  configurations. 

The  exact  solution  for  the  guided-guided  beam  is  given 
as  ir2  EI/L2  for  the  lowest  eigenvalue  P  .  For  nine  internal 
nodes,  the  value  of  X  =  1.0  produced  a  solution  accurate  to 
tt2  within  10  significant  figures.  For  engineering  purposes, 
this  is  the  exact  solution.  When  compared  with  the  convention¬ 
al  solution,  the  trigonometric  finite  differences  were  more 
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Figure  9  demonstrates 


accurate  for  the  range  .71  <  X  <  <*>. 
this  behavior.  The  same  result  occurs  for  19  internal  nodes, 
with  an  improved  accuracy  range  of  .71  <  X  <  “.  Once  again, 
the  optimal  value  of  X  is  1.0.  How  well  does  this  compare  to 
the  mode  shape? 

The  exact  mode  shape  is  W(X/L)  =  costt  (X/L)  .  Once  the 
problem  of  rigid  body  translation  was  accounted  for,  the  eigen¬ 
vectors  for  the  two  mesh  sizes  matched  the  expected  mode  shape. 
This  is  particularly  encouraging  because  no  nodal  displacement 
was  known  ahead  of  time.  In  pinned  and  clamped  beams,  the 
displacement  at  the  boundary  is  known  (equal  to  zero) .  This 
provides  less  unknowns  for  the  program  to  handle.  With  guided 
and  free  boundaries,  the  displacement  is  left  as  an  unknown. 

In  all  guided  and  free  cases  though,  the  computer  returned  the 
correct  displacement. 

From  the  exact  and  generated  mode  shapes,  a  half-buckled 
wavelength  of  L  is  evident.  Therefore,  choosing  X  equal  to 
half  of  the  buckled  wavelength  produces  an  exact  solution. 

General  Results 

Ten  combinations  of  boundary  conditions  were  tested  on  the 
simple  uniform  beam.  In  each  case,  the  exact  buckling  load 
could  be  found  by  ,}^ropeI,,  choice  of  the  wavelength  parameter 
(the  optimal  X,  X  .  Also,  the  eigenvector  exactly  matched 
the  equivalent  points  on  the  known  mode  shape. 

The  number  of  nodes  used  to  model  the  beam  depends  upon 
the  user  and  how  much  time  and  money  are  available.  For  this 
thesis,  time  and  money  had  no  real  effect  on  the  number  of 
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nodes  chosen  because  of  the  simplicity  of  the  beam.  Although 
using  more  nodes  takes  longer,  the  increase  is  not  noticeable. 
The  key  features  to  remember  in  choosing  N  nodes  are  the  in¬ 
crease  in  accuracy  and  the  increase  in  the  size  of  the  gener¬ 
ated  matrices.  Also,  it  is  convenient  to  have  a  node  at  the 
point  of  maximum  deflection.  This  produces  a  smoother  eigen¬ 
vector  shape.  Otherwise,  the  maximum  deflection  will  not  be 
correctly  displayed. 

The  computer  handles  Eq  (4-2)  which  involves  reducing 
matrices  of  order  N  x  N.  Increasing  the  number  of  nodes  will 
increase  the  range  of  X  for  a  good  solution  (compared  to  the 
conventional  solution) ,  improve  the  buckled  mode  shape  by  re¬ 
turning  more  nodal  displacements,  but  will  raise  the  number 
of  elements  in  the  A  and  B  matrices.  For  N  =  9  nodes,  in  a 
pinned-pinned  beam  there  are  81  elements.  With  19  nodes, 
there  are  381  elements.  As  I  stated  before,  the  change  from 
N  =  9  to  N  =  19  has  no  significant  change  on  the  cost  to  the 
user  for  these  simple  beams.  But  the  point  must  be  kept  in 
mind  when  working  with  more  complicated  beams  and  with  complex 
structural  elements  (plates,  shells)  to  which  this  technique 
can  be  applied. 

Nine  nodes  were  chosen  to  produce  the  results  in  Figures 
9  through  18.  Using  nine  nodes  gave  a  convenient  mesh  size 
(.1L)  and  produced  fairly  accurate  buckling  load  plots  and 
eigenvectors.  As  Hannah  demonstrated,  even  a  lower  number  of 
nodes,  such  as  five,  will  produce  an  exact  answer.  However, 
the  choice  of  X  must  be  good  because  the  optimal  range  of  X 
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shrank  by  choosing  fewer  nodes.  So  although  the  exact  answer 
is  available,  the  choice  of  X  can  greatly  affect  your  error. 

For  the  10  beams  analyzed,  the  range  of  X  for  which  Pcr 
was  more  accurate  than  the  conventional  solution  was  checked. 
These  results  are  listed  in  Table  I  under  the  heading  X  range. 
In  all  cases,  the  choice  of  X  which  gave  the  exact  buckling 
load  corresponded  to  half  of  the  buckled  wavelength.  This  is 
the  same  feature  as  seen  in  literature  under  the  effective 
buckling  length  idea  (18). 

TABLE  I 


BOUNDARIES 

*opt 

X 

range 

Pinned-pinned 

1.0 

.71 

< 

X 

< 

oo 

Pinned-clamped 

.632 

.45 

< 

X 

< 

CO 

Pinned-guided 

2.0 

1.42 

< 

X 

< 

oo 

Pinned-free 

1.0 

.71 

< 

X 

< 

oo 

Clamped-clamped 

0.5 

.36 

< 

X 

< 

00 

Clamped-guided 

1.0 

.71 

< 

X 

< 

oo 

Clamped-free 

2.0 

1.42 

< 

X 

< 

00 

Guided-guided 

1.0 

.71 

< 

X 

< 

00 

Guided-free 

2.0 

1.42 

< 

X 

< 

oo 

Free-free 

1.0 

.71 

< 

X 

< 

oo 

If  the  buckling  wavelength  of  the  mode  shape  is  not  known, 

a  choice  of  X  •*  1.4  gives  results  better  than  the  conventional 

approach.  This  result  is  seen  by  comparing  the  ranges  of  X  in 

Table  I,  In  addition,  merely  choosing  a  large  X  will  result  in 

better  accuracy.  For  a  large  X,  the  solution  is  most  likely  a 

lower  bound  on  P  and  therefore  a  conservative  answer. 

cr 

The  general  behavior  of  the  Pcr  vs  X  plots  in  Figures  9 
to  12  explains  the  previous  statement.  It  also  explains  why 
an  exact  solution  is  always  obtainable  for  the  beam  analyzed 
within  the  bounds  of  the  model  used.  Figure  8  typifies  the 
behavior  of  Pcr  vs  X  for  the  uniform  beams. 


Solutions  to  Eq  (4-5)  occur  for 

X  »  0  or  ~  ntr  n  =  0,1,2,...  (4-6) 

The  first  singularity,  at  X  *  0  is  generally  not  a  problem. 
Setting  the  wavelength  equal  to  zero  would  be  of  no  use.  This 
leaves  the  second  part  of  Eq  (4-6) .  The  singularities  generated 
by  it  are 
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(4-7) 


- -r 


3 


A 


(n  =  1,2 


This  forms  a  sequence  of  singularities  from  near  zero  out  to 
A  =  h/2.  For  A  >  h/2,  no  singularities  exist.  This  value  of 
A(h/2)  provides  a  good  start  for  a  Pcr  vs  A  plot.  Below  A  =  h/2, 
singularities  occur.  At  A  =  h/2,  Pcr  =  ».  For  the  uniform 
beam  plot,  the' curves  started  at  A  =  .2  because  the  nearest 
singularity  was  at  A  =  0.05  (h  =  .1).  This  helped  avoid  the 
extremely  large  values  of  the  buckling  load  in  this  region. 

On  the  other  end  of  the  curve  is  the  conventional  finite 
difference  solution.  The  conventional  solution  is  always  a 
lower  bound.  Modeling  the  beam  by  nodes  breaks  the  beam  con¬ 
tinuum  into  a  finite  number  of  degrees  of  freedom,  making  the 
beam  more  flexible  than  the  continuum.  This  increased  flexi¬ 
bility  helps  the  beam  to  buckle  at  a  lower  value  of  Pcr  than 
the  continuum  beam  would.  Therefore,  the  conventional  solu¬ 
tion  is  always  lower  than  the  true  solution. 

Eq  (4-6)  at  n  =  0  implies  A  =  ®  is  a  singularity.  How¬ 
ever,  Eq  (3-29)  showed  that  as  A •*■<»,  h-»-h.  Therefore,  a  very 
large  A  reduces  the  trigonometric  solution  to  the  conventional 
solution  and  not  to  a  singularity. 

The  fact  that  the  exact  solution  is  attainable  is  easily 
explained  now.  At  A  =  h/2,  the  value  of  Pcr  is  positive  in¬ 
finite.  As  A-*»,  the  critical  load  reduces  to  that  determined 
by  the  conventional  solution,  which  is  lower  than  Pcr  exact. 
Therefore,  the  buckling  load  curve  starts  higher  than  the 
exact  solution  and  ends  up  lower  than  the  exact  solution.  By 
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the  mean  value  theorem,  the  curve  must  pass  through  the  exact 
solution  for  some  X,  h/2  <  X  * 

Physically,  increasing  X  results  in  a  buckled  mode  shape 
of  longer  wavelength.  A  smaller  force  is  required  to  buckle 
a  beam  with  this  longer  wavelength.  Therefore,  the  buckling 
load  decreases  with  increasing  X. 

Summary 

By  analyzing  the  simple  uniform  beam,  several  results 
have  been  achieved. 

1)  The  equation  that  generates  the  solution  success¬ 
fully  predicted  the  buckling  load  and  mode  shape 
for  10  known  cases,  for  an  optimal  X. 

2)  The  optimal  X  was  half  of  the  buckled  wave  shape. 

If  the  shape  is  unknown,  a  X  *  1.4  will  give  a 
better  solution  than  the  conventional  finite 
difference. 

3)  The  exact  buckling  load  always  lies  somewhere  on 
Pcr  vs  X  curve. 

4)  The  exact  solution  was  obtained  at  the  same  X  for 
N  -  9  or  19.  However,  N  =  19  provides  a  larger 
range  of  good  X  and  more  points  for  the  eigenvector 
mode  shape . 
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Pig  9.  Buckling  load  curve  comparing  error  for  trig¬ 
onometric  and  conventional  solutions  with  the 
known  solution  (represented  by  the  solid  line) . 
This  plot  is  representative  of  the  pinned-pinned, 
pinned-free,  clamped-guided,  guided-guided,  and 
free-free  beam  solutions.  X  =1.0 
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Fig  10.  Buckling  load  curve  comparing  error  for  trig¬ 
onometric  and  conventional  solutions  with  the 
known  solution  (represented  by  the  solid  line) . 
This  plot  is  representative  of  the  pinned- 
clamped  beam,  ^Qpt  =  *63. 
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Fig  11.  Buckling  load  curve  comparing  error  for  trig¬ 
onometric  and  conventional  solutions  with  the 
known  solution  (represented  by  the  solid  line) . 
This  plot  is  representative  of  the  clamped- 
clamped  beam,  X  =  .5. 
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Fig  12.  Buckling  load  curve  comparing  error  for  trig¬ 
onometric  and  conventional  solutions  with  the 
known  solution  (represented  by  the  solid  line) . 
This  plot  is  representative  of  the  clamped-free 
pinned-guided,  and  guided-free  beams,  A  =  2.0. 
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Fig  13.  Comparison  of  the  calculated  eigenvector  versus 
the  exact  mode  shape.  This  plot  represents  the 
pinned-pinned,  pinned-free,  guided-guided,  and 
free-free  beam. 
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Comparison  of  the  calculated  eigenvector  versus 
the  exact  mode  shape.  This  plot  represents  the 
pinned-guided  or  guided-free  beam. 
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Pig  15.  Comparison  of  the  calculated  eigenvector  versus 
the  exact  mode  shape.  This  plot  represents  the 
clamped-clamped  beam. 
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Comparison  of  the  calculated  eigenvector  versus 
the  exact  mode  shape.  This  plot  represents  the 
pinned-clamped  beam. 
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Fig  17.  Comparison  of  the  calculated  eigenvector  versus 
the  exact  mode  shape.  This  plot  represents  the 
clamped-guided  beam. 
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V.  Beams  With  Variable  Inertia 

Introduction 

The  previous  section  dealt  with  simple  uniform  beams. 

They  were  used  to  verify  the  equation  that  governs  eigenvalue 
determination  as  well  as  the  computer  program  used  to  solve 
that  equation.  The  uniform  beams  tested  had  a  constant  moment 
of  inertia.  If  the  beam's  inertia  can  be  increased  over  its 
middle,  a  higher  load  will  be  required  to  cause  buckling. 

The  beam  is  now  a  more  stable  structure,  compared  to  its 
uniform  cousin.  In  doing  so,  however,  a  beam  with  a  variation 
in  its  moment  of  inertia  results.  Often  the  cross-section 
changes  abruptly,  resulting  in  a  discontinuity  of  the  inertia 
(22).  A  discontinuity  of  this  type  occurs  in  steel  structures, 
where  a  beam  is  strengthened  by  riveting  additional  plates  or 
angles  along  portions  of  the  beam. 

One  of  the  improvements  modeled  in  the  thesis  is  the 
ability  to  handle  changes  in  inertia.  The  moment  of  inertia 
I(x)  is  regarded  as  a  function  of  x,  measured  along  the  beam's 
axis.  The  use  of  I  in  this  manner  includes  it  within  the  in¬ 
tegration  of  the  virtual  work  equation,  Eq  (3-30) .  Inertias 
along  the  beam  are  now  imbedded  within  the  matrix  equation, 

Eq  (3-45) . 

This  section  of  the  thesis  makes  use  of  these  abilities 
by  examining  continuous  and  discontinuous  changes  in  inertia. 
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In  some  cases  the  exact  buckling  load  is  known  beforehand, 
in  others  it  is  only  bounded.  In  each  case,  the  inertia  at 
some  point  along  the  beam  is  defined  as  unity  and  the  remain¬ 
ing  inertias  are  expressed  as  the  ratio  with  respect  to  that 
datum.  The  solution  that  is  generated  is  a  function  of  the 
chosen  inertia.  As  was  done  last  section,  E  and  L  are  defined 
as  one  for  simplicity. 

The  buckling  load  P  versus  buckling  wavelength  X  is 

L  IT 

plotted  for  several  cases  to  check  the  trigonometric  finite 
difference  accuracy,  when  possible.  In  addition,  the  search 
for  the  X  corresponding  to  the  critical  load  will  be  discussed. 

Optimization  of  X 

For  several  of  the  cases  to  be  examined,  only  a  bound 
on  the  solution  is  available  in  the  literature.  The  wave¬ 
length  parameter  corresponding  to  the  buckling  load  bound 

can  be  determined  from  the  P  vs  X  plot  for  the  beam.  But 

cr  r 

that  Pcr  is  not  the  best  solution.  Eventually,  the  technique 
of  this  thesis  will  be  applied  to  buckling  problems  with  no 
known  solution.  What  will  be  the  buckling  load  then? 

In  the  uniform  beam  analysis,  a  range  of  1.4  to  infinity 
for  the  buckling  wavelength  gave  fairly  good  results  for  all 
boundary  conditions.  This  technique  of  looking  for  a  trend 
can  be  extended  by  comparison  with  known  solutions  to  variable 
inertia  beams.  Eventually,  a  rule  of  thumb  for  X  can  be  de¬ 
veloped.  The  corresponding  critical  load  may  be  fairly  good. 

A  promising  procedure  for  determining  an  optimal  buckling 


53 


wavelength  is  put  forward  by  Stein  and  Housner  in  regards  to 
plate  buckling  (2) .  The  technique  is  applicable  to  beam  buck¬ 
ling. 

A  trial  solution  is  performed  for  a  beam  modeled  by  N 
nodes,  equally  spaced  throughout  the  beam.  The  value  of  the 
critical  buckling  load  is  determined  at  several  wavelengths, 
establishing  the  familiar  buckling  load  curve,  P  vs  X. 

The  process  is  repeated  with  a  different  number  of  internal 
nodes  M.  The  buckling  load  curve  is  superimposed  on  the  one 
obtained  for  N  nodes  and  is  shown  in  Figure  18. 


Fig  18.  Optimization  Curves  for  Buckling  Wavelength 
In  the  previous  section,  the  optimal  wavelength  corre¬ 
sponded  to  half  of  the  wavelength  of  buckled  mode  shape.  For 
a  beam  which  is  non-uniform  in  inertia,  the  wavelength  ceases 
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to  have  any  physical  significance.  To  find  the  optimal  wave¬ 
length,  Figure  18  is  made  use  of.  The  value  of  A  at  which 
the  two  curves  cross  indicates  that  convergence  of  the  buck¬ 
ling  load  solution  is  relatively  best  at  this  value  since  for 
this  wavelength,  increasing  the  number  of  nodes  produces  no 
change  in  the  buckling  load.  Increasing  the  number  of  nodes 
will  return  an  improved  solution  at  every  wavelength,  unless 
the  wavelength  corresponds  to  the  exact  buckling  load.  At 

A  in  Figure  18,  increasing  the  number  of  nodes  had  no 
opt 

effect. 

With  this  optimization  method  in  mind,  two  sets  of  beams 
with  a  uniform  mesh,  continuous  and  discontinuous  inertia, 
are  examined. 

Continuous  Inertia 

For  a  beam  with  continuously  varying  inertia,  each  node 
along  the  beam  requires  a  specific  inertia.  Two  examples 
are  discussed  in  which  the  moment  of  inertia  is  defined  by 
a  simple  formula. 

A.  Sinusoidal  Variation.  Brogan  and  Almroth  (18)  give 
a  buckling  solution  for  a  pinned-pinned  beam  with  a  sinusoidal 
variation  in  inertia.  The  beam  illustrated  in  Figure  19, 
has  an  inertia  defined  as 

I(x)  =  I#(l  +  sin1^)  (5-1) 

A  one  term  approximation  to  the  solution  yields  Pcr  = 

1.825ir2EI  /L2.  This  is  higher  than  the  solution  for  a  pinned- 
0 

pinned  beam  with  constant  inertia  I  ,  which  was  found  to  be 
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P__  =  it2EI  /L2  in  the  previous  section,  as  is  expected, 
cr  o 


Fig  19.  Sinusoidal  Variation  in  Inertia 


To  model  the  beam,  nine  internal  nodes  were  used.  The 
x-coordinate  of  each  node  was  calculated,  allowing  use  of 
the  inertia  formula  to  calculate  I (x) .  Using  the  approximate 
solution  as  a  guide,  the  range  .75  <  X  <  °°  provides  answers 
more  accurate  than  the  conventional  finite  difference  solution. 
Because  the  solution  is  not  precisely  known,  the  optimization 
technique  is  tried.  The  first  curve  generated  was  for  nine 
nodes.  The  number  nine  was  chosen  because  it  places  a  node 
at  the  position  of  maximum  inertia,  in  much  the  same  way  as 
a  node  at  maximum  deflection  of  the  uniform  beam  was  needed 
to  fully  model  the  mode  shape.  With  this  concept  as  a  guide 
for  choosing  another  set  of  nodes,  the  choice  of  N  =  11  was 
made.  Any  odd  number  of  internal  nodes  will  place  a  node  at 
the  center.  Although  this  requires  a  second  solution  to  the 
eigenvalue  equation,  the  size  of  the  matrix  is  not  much  larger 
than  for  N  =  9.  Size  depends,  in  general,  upon  N2  as  was 
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shown  earlier.  The  two  matrix  sizes  combined  are  still  small¬ 
er  than  one  corresponding  to  N  =  19  and  may  produce  a  better 
solution. 

For  N  =  11,  the  range  of  wavelengths  for  which  the  trig¬ 
onometric  series  is  more  accurate  shrinks  to  .8  <  X  <  <». 

The  shrinkage  is  expected.  As  N  increases,  the  conventional 
solution  improves  slightly  faster  than  does  the  trigonometric 
approach.  However,  to  get  the  conventional  solution  very 
close  to  the  buckling  load  requires  many  nodes,  whereas  the 
trigonometric  solution  requires  only  a  good  choice  of  wave¬ 
length  with  fewer  nodes. 

The  curves  generated  for  N  =  9  and  11  are  shown  in 

Figure  23.  They  intersect  at  X  =  1.018.  The  corresponding 

buckling  load  is  1.8268ir2EI  /L2 ,  slightly  higher  than 

o 

lower  bound  using  Rayleigh  -  Ritz.  For  this  non-uniform 

pinned-pinned  beam,  the  optimization  technique  worked  well. 

Two  more  solutions  were  obtained  using  10  and  12  nodes. 

These  nodal  arrangements  do  not  model  the  maximum  moment  of 

inertia.  Directly  comparing  the  buckling  curves  produced 

for  N  =  10  and  12  provides  a  buckling  solution  of  1. 827it2EI#/L2 

at  X  =  1.01.  Comparing  the  curves  for  N  =  9  and  10,  which 

have  on-  and  off-center  nodes,  results  in  P  =  1.8267tt2EI  /L2 

cr  o 

at  X  =  1.02.  These  results  show  that  an  accurate  solution  can 
be  found  if  both  curves  use  the  same  nodal  arrangement,  i.e. 
both  curves  have  on-center  or  off-center  nodes,  in  this  example. 
Comparing  dissimilar  nodal  arrangements  (N  *  9,  10)  resulted 
in  a  fairly  good  solution  also.  This  was  not  expected  and  will 
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be  checked  with  the  next  model. 

An  additional  nodal  arrangement,  N  =  5,  was  checked. 

The  range  of  improved  accuracy  over  the  conventional  solution 
increased  to  .73  <  X  <  °°.  Comparing  the  curve  of  Pcr  vs  X 
with  another  on-center  solution  (N  =  9)  gives  the  same  solu¬ 
tion  as  with  N  =  9  and  11.  Therefore,  in  some  cases,  fewer 
nodes  could  be  used  to  predict  the  critical  load. 

One  final  note  on  the  sinusoidal  pinned-pinned  beam. 

As  discussed  in  the  introduction  to  this  section,  a  beam 
built-up  in  the  center  would  buckle  at  a  higher  load  than 
its  uniform  counterpart.  This  is  seen  for  the  sinusoidal 

beam  which  buckles  at  about  1.827tt2EI  /L2  versus  1.0ir2EI  /L2 

0  0 

for  its  uniform  inertia  cousin.  What  is  surprising  though 
is  that  the  ranges  for  a  good  choice  of  X  and  the  optimal 
wavelength  itself  are  close  to  the  uniform  pinned-pinned  beam. 
Perhaps  a  uniform  beam  can  serve  as  a  good  guide  in  the  choice 
of  X.  This  would  be  advantageous  because  results  for  all  10 
combinations  of  boundary  conditions  are  well  known. 

B.  Exponential  Variation.  Sinha  and  Chou  (17)  present 
several  solution  methods  for  the  axial  buckling  load  of  a 
pinned-pinned  beam  with  an  exponentially  varying  moment  of 
inertia.  The  beam,  similar  to  that  illustrated  in  Figure  19, 


has  an  inertia  defined  as 

I(x)  =  I  e2x/L 
o 


for  0  _<  x  <  — 


(5-2) 


The  beam  is  symmetric  about  the  mid-section. 

The  exact  solution  for  the  axial  buckling  load,  using 
the  same  assumptions  embodied  in  this  thesis,  is  given  in 
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terms  of  the  Bessel  functions. 


J  (n)Y  (ne-J*)  -  J  (ne-J*)Y  (n)  =  0  (5-3) 

0  1  10 

where 

n2  =  P  L2/EI  (5-4) 

cr  o 

Solution  of  Eq  (5-3)  yields  the  beam  buckling  load  as 

1.9634ir2EI  /L2 .  This  is  a  higher  buckling  load  than  was 
o 

found  for  the  sinusoidal  variance  because  the  inertia  is 
larger  in  the  beam  center. 

The  beam  is  modeled  with  an  odd  number  of  nodes  so  that 
the  maximum  inertia,  at  the  beam's  center,  is  included.  The 
first  mesh  tried  was  N  =  9.  The  inertia  at  each  node  was 
calculated  using  Eq  (5-2) .  The  wavelength  range  for  which 
the  trigonometric  solution  is  more  accurate  than  the  conven¬ 
tional  solution  is  1.04  <  X  <  <*>,  with  the  calculated  solution 
at  X  =  1.475. 

Although  the  exact  solution  is  known  and  helped  verify 

the  ability  to  generate  an  accurate  critical  load,  it  provides 

an  opportunity  to  try  the  two  curve  optimization  again.  A 

second  nodal  arrangement,  N  =  11,  is  applied  to  the  problem. 

When  superimposing  the  N  =  9  and  11  curves  in  Figure  24,  the 

intersection  occurs  at  X  =  1.48  with  a  corresponding  solution 

of  P  =  1.9628tt2EI  /L2.  This  solution  is  within  .031%  of  the 
cr  o 

calculated  solution  of  Eq  (5-3)  and  is  most  likely  a  combina¬ 
tion  of  error  in  solving  the  eigenvalue  problem  and  solving 
for  n  in  Eq  (5-4).  A  change  in  n  of  -.014%  would  give  the 
solution  returned  by  the  intersection  of  the  two  curves.  Be¬ 
cause  the  tables  used  to  solve  Eq  (5-3)  were  accurate  for  n 
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only  to  one  part  in  a  thousand,  the  exact  solution  cannot  be 
calculated  as  closely  as  the  intersection  did.  Therefore, 
the  two  curve  intersection  gives  what  can  be  termed  an  exact 
solution. 

As  was  done  for  the  sinusoidal  beam,  two  solution  curves 
corresponding  to  N  =  10,  12  (off-center  nodes)  were  generated. 

The  intersection  of  these  two  curves  produced  P_  =  1.963ir2EI  /L2 

cr  o 

at  X  =  .88,  which  is  fairly  close  to  the  result  from  N  =  9,  11 
except  the  buckling  wavelength  is  quite  different.  The  major 
difference  from  the  sinusoidal  beam  is  the  lack  of  an  inter¬ 
section  between  the  curves  for  N  =  9,  10  and  N  =  9,  12.  These 
curves  are  shown  in  Figures  25  and  26.  The  sinusoidal  beam 
did  return  a  solution  for  the  buckling  load  when  comparing 
meshes  with  nodes  on-  and  off-center.  The  exponential  beam 
does  not  produce  this  result.  Normally,  conventional  solution 
accuracy  improves  as  more  nodes  are  used  because  the  beam  con¬ 
tinuum  is  modeled  by  more  points.  Therefore,  increasing  the 
number  of  nodes  will  raise  the  conventional  solution,  which  is 
the  asymptote  of  the  buckling  load  curve.  If  the  number  of 
nodes  increases,  and  the  conventional  solution  improves,  the 
two  buckling  curves  will  intersect.  Failure  to  model  the  same 
points  of  interest  may  result  in  a  decrease  in  the  conventional 
solution  for  an  increase  in  the  number  of  nodes.  The  buckling 
load  curves  will  not  intersect  in  this  case.  For  nine  nodes, 
the  trigonometric  buckling  load  curve  approached  Pcr=  1.955ir2, 
the  conventional  solution.  For  N  =  10  and  12,  the  buckling  loads 
are  1.946tt2  and  1.951ir2  respectively.  Comparing  N  =  9  and  10 
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or  N  =  9,  12  shows  the  decrease  in  load  for  an  increase  in 
nodes.  Therefore,  no  intersection  will  result.  There  is  an 
increase  in  the  conventional  load  from  N  =  10  to  12,  account¬ 
ing  for  the  solution  mentioned  earlier  for  the  two  meshes. 

An  intersection  with  N  —  9  may  be  possible  for  N  =  14,  if  the 
conventional  solution  is  greater  than  for  N  =  9.  Looking  at 
the  increase  from  N  =  10  to  12  indicates  just  such  a  possibil¬ 
ity.  However,  the  number  of  elements  in  the  eigenvalue  matrix 
is  increasing  rapidly,  bringing  up  the  cost  of  a  solution.  It 
would  be  better  just  to  compare  meshes  which  model  the  same 
points  of  interest. 

For  the  exponential  beam,  inertias  ranged  from  I  to 

0 

2.721  ,  making  it  stiff er  than  its  uniform  cousin  with  inertia 
o 

I  .  According,  the  exponential  beam  buckled  at  a  higher  load. 
0 

The  buckling  wavelength  could  be  greater  or  less  than  the 
uniform  beam  solution  (X  =  1.0)  depending  upon  the  mesh  used. 
This  points  out  the  advantage  of  using  the  optimization  tech¬ 
nique  for  finding  the  buckling  load.  When  a  mesh  was  used 
which  modeled  the  maximum  inertia,  a  buckling  wavelength  was 
obtained  which  was  higher  than  for  the  less  stiff  uniform 
beam. 

Discontinuous  Inertia 

A  beam  with  increased  inertia  about  its  center  provides 
greater  resistance  to  buckling.  To  investigate  variations  in 
inertia,  continuously  changing  beams  were  examined.  But  the 
change  in  cross-section  is  usually  abrupt,  mostly  due  to 
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welding  or  riveting  material  at  strategic  locations  along  the 
column.  The  abrupt  change  leads  to  an  investigation  of  dis¬ 
continuous  inertias. 

Figure  20  illustrates  one  type  of  discontinuous  beam. 


P  =  mEI_/L2 
cr  2 


Fig  20.  Typical  Discontinuous  Beam 

The  beam  is  clamped  at  one  end  and  is  free  at  the  other.  It 

is  discontinuous  in  inertia  at  x  =  aL.  The  buckling  load  is 

defined  in  terms  of  the  modulus  E,  the  total  beam  length 

L,  and  the  inertia  1^.  For  simplicity,  1^  is  defined  as  one. 

The  inertia  1^  will  be  expressed  as  the  ratio  of  I^/I^  with 

I  =1.0.  The  non-dimensional  coefficient  m  is  well  tabulated 
2 

in  literature  (7,8,22,24).  Because  this  type  of  discontinuous 

beam  has  known  solutions  widely  available,  it  will  be  the 

model  for  my  discontinuous  inertia. 

As  was  done  with  the  continuously  varying  beams,  the 

inertia  must  be  calculated  and  provided  for  each  node.  This 

is  not  too  difficult,  there  being  only  two  values,  I  and  I  , 

1  2 

in  the  example  considered.  Past  experience  has  shown  that 
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significant  positions  along  the  beam  should  be  modeled  with 
a  node.  These  have  included  maximum  inertia  and  maximum  dis¬ 
placement.  With  this  in  mind,  a  node  will  be  positioned  at 
x  «  aL,  the  discontinuity  of  the  beam.  Positioning  the  node 
at  that  point  forces  certain  mesh  sizes  on  the  modeling.  For 
the  beams  of  interest,  a  is  a  multiple  of  0.1.  Node  numbers 
of  9,  14,  and  19  are  useful  in  modeling  a  discontinuity  at 
this  position.  The  effect  of  a  node  not  at  the  discontinuity 
will  also  be  examined. 

If  a  node  is  located  at  a  discontinuity,  what  value  of 
inertia  will  be  assigned  to  it?  Should  the  inertia  to  its 
left  (IL)  or  to  its  right  (IR)  be  used?  Possibly  a  simple 
average  or  a  weighted  average  should  be  used. 

Girijavallabhan  (14)  suggests  an  effective  inertia  for 
a  discontinuity  at  node  i 


t  2  + 


Wi-1 


1  Wifl 

+  !R  2  +  wi 


(5-5) 


for  the  beam  of  Figure  20.  This  solution  makes  use  of  an 
assumed  mode  shape  which  is  not  available  in  the  approach  of 
this  thesis.  An  adequate  value  for  I  could  be  found  by  assuming 
that  the  displacements  are  the  same  as  for  the  uniform  beam. 

This  inertia  could  be  used  to  calculate  a  better  eigenvector 
which  in  turn  would  generate  a  better  effective  inertia.  The 
cycle  would  continue  until  some  type  of  convergence  is  obtained. 
The  final  value  for  I  would  be  used  at  the  discontinuity. 

A  simpler  solution  is  provided  by  O'Rourke  (24).  The 
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effective  moment  of  inertia  at  a  step  change  is  given  by 


2IHTL 

eff  “  (lp  +  IT) 


(5-6) 


LR  '  ~L' 

O’Rourke  states  that  this  value  of  the  inertia  represents 
a  value  between  the  inertia  for  a  pinned-pinned  beam  using 
Rayleigh's  method  with  a  step  change  and  the  inertia  for  a 
clamped  beam.  The  real  value  of  Eq  (5-6)  was  not  apparent  in 
his  paper.  However,  Ghali  and  Neville  (7)  derive  this  ex¬ 
pression,  which  is  repeated  in  Appendix  C.  In  doing  so,  con¬ 
tinuity  in  slope  and  moment  is  maintained,  two  very  important 
properties.  This  expression  maintains  the  necessary  proper¬ 
ties  which  is  why  Eq  (5-6)  is  used  in  this  thesis. 


Specific  Examples 

As  a  check  on  the  ability  to  predict  the  correct  eigen¬ 
value  (buckling  load) ,  many  beams  of  the  type  shown  in  Figure 
20  were  tested.'*  The  parameter  "a"  was  varied  from  0.1  to  0.9 

with  a  similar  range  in  I  /I  .  In  all  cases,  using  the  inertia 

1  2 

of  Eq  (5-6)  at  the  discontinuity,  the  eigenvalue  equation  gen¬ 
erated  the  known  solution  for  a  specific  buckling  wavelength. 

As  an  example,  the  following  values  were  chosen  for  the  beam 
shown  in  Figure  21: 

a  ■  .7L 

I  /I  =  .3 
1  2 

The  exact  solution,  as  shown  in  Figure  21,  is  P__  =  2.221016EI 

2 

/L 2 .  The  inertia  calculated  to  represent  the  discontinuity  is 
.4615  using  Eq  (5-6).  For  this  clamped-free  beam,  a  range  of 
. 78  <  A  <  »  provides  accuracy  better  than  the  conventional 


[•Til 
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Fig  21.  Discontinuous  Beam  of  the  First  Example 


A  second  example  looks  more  closely  at  the  effect  of  the 

number  of  nodes.  Choosing  the  values  of  a  =  .4  and  I  /I  =  .4 

1  2 

gives  an  exact  solution  of  M  =  1.66931.  A  Pcr  vs  X  curve  for 
N  =  9  is  generated  with  an  accuracy  range  of  .99  <  X  <  <»  and 
an  optimal  solution  at  X  =  1.4.  A  second  curve,  using  N  =  19 
is  superimposed  in  Figure  27.  The  two  curve  intersection 
occurs  at  X  =  1.4  as  expected. 

Two  more  curves  are  generated,  one  for  N  =  8  and  another 
for  N  *  10.  These  nodal  arrangements  do  not  position  a  node 
at  the  discontinuity.  Therefore,  a  node  sees  only  an  inertia 
of  1.0  or  0.4,  depending  upon  its  location.  For  the  eight 
internal  node  beam,  X  =  .78  produces  the  known  solution.  For 
10  nodes,  though,  the  known  solution  is  never  reached.  This 
buckling  load  curve  appears  in  Figure  28.  The  conventional 
solution,  which  is  the  limit,  is  greater  than  the  known  solu¬ 
tion.  This  could  be  due  to  the  size  of  the  inertia  gradient 
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across  the  discontinuity.  Assuming  a  linear  gradient  across 
the  discontinuity  results  in  effective  inertias  of  .64  and 
.76  for  the  8  and  10  node  meshes  respectively.  The  effective 
inertia  from  Eq  (5-6)  is  approximately  .57.  While  this  is 
the  value  which  maintains  continuity,  there  may  be  an  upper 
limit  to  the  effective  inertia  below  which  the  correct  buck¬ 
ling  load  is  returned  for  some  value  of  the  wavelength.  In 
this  case,  the  10  node  mesh  models  a  beam  that  is  too  stiff. 

As  a  check  on  where  the  upper  limit  on  effective  inertia 
lies,  several  values  for  the  inertia  at  the  discontinuity 
were  substituted.  The  indicator  that  an  effective  inertia  is 
too  high  is  a  conventional  finite  difference  solution  that  is 
greater  than  the  known  solution.  The  beam  is  then  too  stiff. 

An  average  of  IL  and  IR,  giving  the  effective  inertia 
as  .7,  was  tested.  For  all  values  of  the  buckling  wavelength 
the  known  solution  was  not  produced.  Therefore,  the  10  node 
inertia  of  .76,  being  higher  than  .7,  would  be  expected  to 
be  ineffective. 

A  weighted  average  was  also  used  to  produce  an  effective 
inertia.  This  procedure  sums  the  inertia  of  each  segment  - 
weighted  by  the  percentage  of  the  beam  it  represents  -  and 
divides  by  the  sum  of  the  inertias.  The  resulting  inertia  of 
.457  is  below  the  value  returned  by  Eq  (5-6),  .5714.  The 
weighted  average  produced  a  correct  solution  at  A  =  .44. 

While  this  method  produced  a  larger  range  of  wavelengths  for 
which  the  trigonometric  solution  was  better  than  the  conven¬ 
tional  solution  than  for  the  inertia  of  Eq  (5-6) ,  the  real 
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results  were  not  as  good.  The  conventional  solution  is  worse 
for  a  weighted  inertia  because  continuity  is  not  maintained. 
So  while  you  may  have  a  larger  range  of  wavelength  at  which 
the  solution  is  better  than  conventionally,  the  conventional 
solution  itself  is  not  very  good.  A  better  test  is  how  large 
a  bandwidth  on  X  will  give  an  accuracy  within  ±.1%  of  the 
known  buckling-  load.  For  a  weighted  inertia,  this  bandwidth 
extends  from  wavelengths  equal  to  .43  to  .48,  a  fairly  small 
bandwidth.  The  inertia  obtained  using  continuity  gives  a 
bandwidth  from  1.26  to  1.61,  five  times  larger.  In  fact,  an 
inertia  from  Eq  (5-6)  gave  a  larger  range  of  bandwidth  accu¬ 
racy  than  for  any  other  inertia.  In  addition,  there  is  an 
upper  limit  on  the  effective  inertia  above  which  the  correct 
solution  is  unobtainable.  This  point  should  be  kept  in  mind 
when  modeling  without  a  node  at  a  discontinuity. 

In  addition,  comparing  curves  for  N  =  8  and  N  =  9  yields 

an  intersection  for  a  buckling  load  of  1.75EI  /Lz,  far  above 

2 

the  known  solution  of  1.669;  in  Figure  29.  This  demonstrates 
the  potentially  poor  results  of  trying  to  optimize  the  wave¬ 
length  with  meshes  which  do  not  model  the  same  points. 

A  study  should  be  made  of  what  is  the  best  effective 
inertia  and  compare  it  to  that  derived  by  Eq  (5-6).  For  the 
purposes  of  this  thesis,  Eq  (5-6)  has  been  shown  effective 
for  all  discontinuity  combinations  tested. 

A  final  example  of  a  more  complex  discontinuous  beam  is 


( 


shown  in  figure  22.  A  pinned-pinned  beam  with  two  discontin¬ 
uities  is  examined.  The  effective  inertias  at  x  =  ,2L  and 


x  *  .8L  are  determined  by  Eq  (5-6).  Using  nine  internal  nodes, 
so  chosen  to  place  a  node  at  each  discontinuity,  provides  the 
expected  solution  at  X  =  .64.  A  second  nodal  arrangement, 

N  =  19,  is  used  to  intersect  with  the  curve  produced  by  N  =  9. 


Fig  22.  Third  Example  of  a  Discontinuous  Beam 
The  intersection  occurs  for  X  =  .64,  P  =  8.51,  matching  the 

Ci 

known  solution.  The  results  are  shown  in  Figure  30. 

Summary 

A  beam  with  a  discontinuous  or  varying  moment  of  inertia 
can  be  analyzed  using  the  approach  of  the  thesis.  In  all  such 
beams,  an  attempt  should  be  made  to  place  a  node  at  the  maximum 
moment  of  inertia  and  also  at  any  discontinuities.  For  a  con¬ 
tinuous  variation  in  inertia,  the  uniform  beam  results  can  pro¬ 
vide  an  idea  of  a  good  choice  of  X,  generally  larger  than  the 
uniform  solution.  The  discontinuous  beams  can  also  use  this 


principle. 

An  accurate  method  of  determining  the  correct  value  for 


the  buckling  load  is  to  plot  P  vs  X  for  two  sets  of  nodes. 
The  intersection  of  these  two  curves  is  the  desired  solution 
The  two  sets  of  nodes  should  be  somewhat  similar.  Also,  the 
effective  inertia  at  a  discontinuity  is  evaluated  using  Eq 
(5-6)  . 

With  these  principles,  the  known  solution  was  derived 
for  every  example  tested. 
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Fig  24,  Buckling  load  curves  for  a  beam  with 
exponential  inertia.  N  =  9*11. 
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Fig  26.  Buckling  load  curves  for  a  beam  with  exponen 
tial  inertia.  N  =  9,12. 


P-CRIT 


Fig  27.  Buckling  load  curves  for  a  discontinuous 
beam,  lx/lz  =  .4,  a/L  =  .4.  N  =  9,19 

P  *  1.669  EI./L2,  \  *  1.4. 
cr 
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Fig  28.  Buckling  load  curve  for  a  discontinuous  beam, 

Ij/I2  =  .4,  a/L  =  .4.  N  =  8,9. 

P  *  1.669  EI./L2,  X  =  1.4 
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Fig  29.  Buckling  load  curves  for  a  discontinuous, 
Ix/I2  *  .4,  a/L  =  .4.  N  =  8,9. 

Pcr  =  1.669  EI2/L2,  A  =  1.4 
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Fig  30.  Buckling  load  curves  for  a  discontinuous  beam 
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N  *  9,19.  P  *  8.51,  X  =  .64 
cr 


VI.  Variable  Mesh 

Introduction 

In  the  previous  section,  beams  with  various  inertias 
were  examined.  Often,  the  beam  had  a  geometric  discontinuity 
in  its  inertia.  This  discontinuity  was  modeled  with  a  node 
at  the  location  of  the  step  change  to  produce  the  best  results 


Increasing  the  number  of  nodes  also  improved  the  solution 
accuracy,  if  the  discontinuity  is  modeled.  Previous  examples 
showed  the  increase  in  accuracy  between  nine  and  nineteen 
nodes,  both  of  which  placed  anode  at  a  discontinuity.  The 
change  to  19  nodes  decreased  the  mesh  size  by  a  factor  of  two, 
increasing  the  solution  accuracy. 

To  obtain  increased  accuracy  up  until  now,  the  number 
of  nodes  had  to  be  increased.  Doing  this  increases  the  time 
and  cost  to  the  user.  An  alternative  approach  to  increasing 
the  number  of  nodes  is  to  respace  the  available  nodes.  For 
the  discontinuous  beams,  a  nine  node  solution  may  give  better 
accuracy  if  the  available  nodes  are  respaced.  With  several 
nodes  at  and  around  the  discontinuity,  a  better  range  of 
critical  buckling  loads  may  be  obtained  than  is  possible  with 
a  uniform  mesh  of  the  same  number  of  nodes.  Figure  31  shows 
the  difference  in  node  spacing  for  a  discontinuous  beam  with 
nine  nodes  for  uniform  and  variable  meshes. 
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Fig  31.  Comparison  of  Variable 
Uniform  and  Variable  Meshes 

The  variable  mesh  has  taken  the  nine  nodes  with  uniform 
spacing  and  redistributed  several  of  them  about  the  discon¬ 
tinuity.  The  remaining  nodes  are  spaced  according  to  the 
sizes  of  the  remaining  beam  sections.  Placing  more  nodes 
nearer  the  discontinuity  than  is  available  from  a  uniform 
mesh  hopefully  gives  the  same  effect  as  increasing  the  total 
nodes  used.  For  example,  a  nine  node  beam  has  an  h  =  .1  mesh 
size.  If  two  nodes  are  moved  closer  to  the  discontinuity  for 


a  mesh  spacing  of  h  =  .05  (Figure  31),  the  mesh  size  near  the 

2 

change  is  the  same  as  for  N  =  19.  The  question  that  must 
be  answered  is  how  the  remaining  nodes  can  be  distributed 
for  a  more  accurate  solution  to  result.  If  this  is  possible, 
better  accuracy  can  be  obtained  by  clustering  more  nodes  about 
discontinuities,  points  of  maximum/minimum  inertia,  and  other 
significant  locations.  Section  VI  examines  just  this  question. 

Technique 

Section  III  of  this  thesis  showed  the  numerical  integra¬ 
tion  technique  used  to  solve  the  virtual  work  expression  with 
a  uniform  mesh.  With  a  variable  mesh,  the  numerical  integra¬ 
tion  technique  must  be  expanded  to  accomodate  changes  in  mesh 
size  throughout  the  beam.  Figure  32  shows  this  scheme  for  a 
beam  with  three  different  mesh  sizes.  This  figure  is  the 
analog  to  Figure  5  which  works  for  uniform  mesh.  Nodes  i  and 
j  correspond  to  the  nodes  at  which  mesh  size  changes.  The 
three  mesh  sizes  are  h  ,  h  ,  h  .  With  the  change  in  mesh 

12  3 

size,  Eq  (3-33)  must  be  reintegrated.  Applying  the  trape¬ 
zoidal  rule,  the  equation 

L  L 

/  fdx  =  P  /  gdx  (3-33) 

0  0 
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in  Eq  (6-1).  Each  bracketed  ({})  set  of  terms  has  the  same 
construction  as  the  others,  except  for  the  different  mesh 
size  and  end  points.  It  becomes  fairly  simple  to  program  this 
repetitious  pattern.  The  end  result  is  the  construction  of 
the  eigenvalue  matrices  for  each  segment  of  uniform  mesh  size. 
These  eigenvalue  matrices  correspond  to  the  ones  described  by 
equation  (3-45) .  However,  whereas  the  matrices  generated  by 
Eq  (3-46)  were  global  -  describing  the  entire  beam,  the  matri¬ 
ces  built  by  a  segment  of  constant  mesh  size  h^  are  'local' 

-  corresponding  to  only  that  segment.  To  produce  the  required 
global  matrix  equation,  the  local  matrices  must  be  added  for 
each  segment  of  constant  mesh  size.  For  the  beam  in  Figure 
32  ,  there  are  three  segments  of  constant  mesh  size.  There¬ 
fore,  three  local  eigenvalue  matrices  are  produced  correspond¬ 
ing  to  the  bracketed  terms  in  Eq  (6-1) .  These  matrices  are 
added  together,  making  sure  the  nodes  overlap  correctly,  pro¬ 
ducing  a  global  matrix  equation  of  the  same  form  as  Eq  (3-46) . 
This  process  is  similar  to  the  construction  of  a  global  stiff¬ 
ness  matrix  in  finite  elements,  in  which  segments  of  a  beam 
add  their  local  stiffness  matrices  together  to  form  a  global 
matrix.  The  local  matrices  developed  here  depend  upon  stiff¬ 
ness,  inertia,  and  the  displacements  at  each  node,  as  well 
as  the  mesh  size. 

Because  various  mesh  sizes  are  used  in  the  beam,  diffi¬ 
culties  arise.  Eq  (6-1)  demonstrates  the  need  to  define 
f^_j^,  the  combination  of  derivatives  at  a  node  where  the  mesh 
changes. 


82 


Previously, 
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This  expression  used  a  full-station  central  difference  which 
was  derived  using  a  constant  mesh  size.  But  Figure  32  shows 

that  f.  .  has  mesh  sizes  h  and  h  on  either  side  of  node  (i-1) . 

1-1  i  2 

The  central  difference  expressions  created  for  the  second  de¬ 
rivative  in  Section  III  are  of  no  use  at  a  node  separating 
two  different  mesh  sizes. 

Using  a  technique  similar  to  the  one  used  in  Section  III, 
a  finite  difference  expression  for  the  second  derivative  is 
calculated  in  Appendix  B.  At  a  mesh  discontinuity  node,  this 
expression  is  used  instead  of  the  second  derivative  shown  in 
Eq  (6-2) .  The  first  order  derivatives  do  not  have  this  problem 
because  their  reference  points  are  between  two  nodes,  where 
the  mesh  does  not  change  size.  With  the  appropriate  nodal 
expressions  used,  the  eigenvalue  matrix  is  calculated  in  the 
same  manner  as  in  Eq  (3-46). 

The  mesh  discontinuities  introduce  one  more  unfavorable 
aspect  to  the  problem.  Previously,  in  Eq  (3-34) ,  the  mesh 
size  h  was  factored  out  of  the  eigenvalue  equation  because  it 
was  constant  and  appeared  on  both  sides  of  the  equation.  With 
the  variable  mesh,  the  mesh  size  for  each  segment  Ik  must  be 
retained  in  the  problem.  Of  greater  significance  is  the  need 
to  multiply  each  term  of  the  A  matrix  of  Eq  (3-46)  by  1/h? 
where  Ik  is  the  equivalent  mesh  corresponding  to  each  segment 
of  the  beam.  For  a  constant  mesh  size,  the  equivalent  mesh  h 
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could  be  factored  outside  the  A  and  B  matrices  in  Eq  (4-1) 
and  (4-2).  Doing  so  permitted  solving  the  eigenvalue  matrix 


for  k,  Eq  (4-3),  then  varying  X  to  obtain  P.,..  Because  h  is 
not  constant  throughout  the  beam,  neither  is  h.  The  equivalent 
mesh  term  is  now  brought  within  the  matrices  of  Eq  (3-46) ,  and 
with  it  goes  the  ability  to  vary  X  after  the  eigenvalue  matrix 
problem  is  solved.  The  equivalent  mesh  sizes  cannot  be  ex¬ 
pressed  as  ratios  of  one  another,  in  an  effort  to  retain  the 
variation  in  X  external  to  the  matrix  solution,  because  the 
ratios  themselves  depend  upon  X.  The  mesh  sizes  h^  could  be 
scaled  this  way,  but  that  is  of  no  significant  use. 

The  use  of  a  variable  mesh  requires  that  the  wavelength 
parameter  be  specified  before  solving  the  eigenvalue  equation. 
The  eigenvalue  returned  is  the  buckling  load  and  not  k  in 
Eq  (4-3) ,  which  simplified  solving  problems  with  a  uniform 
mesh.  Because  of  this,  a  Pcr  vs  X  curve  cannot  be  generated 
with  only  one  solution  of  the  matrix  equation.  The  variable 
mesh  approach  requires  constructing  and  solving  the  matrix 
equation  each  time  a  point  on  the  buckling  load  curve  is  de¬ 
sired.  Before,  the  buckling  load  curves  in  Sections  IV  and  V 
required  one  matrix  solution.  Now,  twenty  solutions  may  be 
needed  to  produce  the  same  curve.  The  increase  in  time  and 
cost  makes  the  variable  mesh  approach  appear  unusable  right 
away.  Only  for  a  better  range  of  accuracy  than  is  available 
with  an  equal  number  of  nodes  in  a  uniform  mesh  can  the  tech¬ 
nique  hope  to  be  effective. 


Test  Results 

To  answer  the  question  just  posed,  whether  variable 
meshes  are  more  accurate  than  uniform  meshes  for  an  equal 
number  of  nodes,  beams  from  Section  V  were  reanalyzed.  These 
beams  have  a  varying  moment  of  inertia,  providing  the  best 
test  of  preferentially  grouped  nodes. 

The  first  beam  tested  appears  in  Figure  33. 


Fig  33.  First  Variable  Mesh  Example 
The  known  exact  solution  is  P  =  .366875  El  /L2 ,  occurring 
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at  X  =  1.61  for  nine  or  19  nodes  uniformly  spaced.  The  number 
of  nodes  to  be  used  for  the  variable  mesh  is  nine.  Nine  nodes 
provides  a  convenient  mesh  size  (.1L)  and  places  a  node  at  the 
discontinuity.  Because  the  discontinuity  at  x/L  =  .2  is  the 
region  of  significance,  more  nodes  will  be  grouped  around  it. 
There  will  be  two  uses  made  with  the  variable  mesh,  one  with 
three  nodes  near  the  discontinuity  and  one  with  five  nodes. 

The  neshes  will  be  used  to  generate  the  buckling  load  curve 
P__  vs  X  (by  multiple  computer  runs)  which  can  be  compared  to 


the  curve  for  nine  nodes  uniformly  generated 


The  first  variable  mesh  tested,  (mesh  A) ,  consists  of 
three  nodes  clustered  more  closely  around  the  discontinuity 
than  is  available  with  the  uniform  mesh.  The  nodal  arrange¬ 
ment  is  shown  in  Figure  34.  The  mesh  size  around  the  dis¬ 
continuity,  h^,  was  chosen  as  one-half  the  uniform  mesh  size 
of  .  1L.  This  corresponds  to  modeling  the  region  near  the  dis- 
continuity  with  a  19  node  mesh,  which  would  provide  a  mesh  of 
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Fig  34.  Discontinuous  Beam  With  Variable  Mesh  A 
The  six  remaining  nodes  (nine  minus  three)  are  distribut¬ 
ed  by  the  percentage  of  remaining  beam.  This  results  in  mesh 
sizes  of  h  =  .075L  and  h  =  .125L.  With  the  first  variable 

I  3 

mesh  in  place,  the  buckling  load  was  calculated  for  20  values 

of  the  wavelength  parameter  (.1  to  2.0) 

A  second  mesh  (B)  was  built  using  five  nodes  about  the 

discontinuity,  also  with  mesh  size  h  =  . 05L.  The  remaining 

2 

four  nodes  are  distributed  to  create  mesh  sizes  h  ■  .10  and 


=  .14,  shown  in  Figure  35. 


tions  of  the  eigenvalue  equation. 

The  curves  obtained  for  meshes  A  and  B  are  displayed 

with  the  nine  node  uniform  mesh  curve  in  Figure  37.  At 

X  -  1.6,  the  uniform  mesh  provides  a  solution  of  P  = 

. 36688EI  /L2  which  is  in  error  by  only  .0014%  of  the  known 
2 

solution.  The  variable  meshes  A  and  B  on  the  other  hand  have 
solutions  of  .36720  and  .36753,  respectively.  The  errors  for 
the  variable  meshes  are  .089%  and  .179%,  much  larger  than  for 
the  uniform  mesh  solution.  In  fact,  figure  37  shows  that  for 
all  X  <  1.6,  the  uniform  solution  is  more  accurate  than  either 
variable  mesh  solutions.  Likewise,  three  nodes  around  the 
discontinuity  are  better  than  five.  For  the  range  X  >  1.6, 
there  are  small  intervals  in  which  each  variable  mesh  solution 
is  more  accurate  than  the  uniform  solution.  Not  because  of 
the  better  quality  of  the  mesh,  but  because  they  fail  to  match 
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the  uniform  curve.  Mesh  A  arrives  at  the  exact  solution  at 
X  =  1.73  and  mesh  B  at  X  =  1.78.  During  the  interval 
1.7  <  X  <  1.9,  the  two  meshes  are  more  accurate.  But  this 
accuracy  quickly  disappears  and  soon  the  uniform  solution  is 
best  again. 

As  the  wavelength  parameter  grows  very  large,  the  uniform 
solution  approaches  the  conventional  solution.  Unfortunately, 
the  variable  mesh  solutions  pass  right  through  this  uniform 
mesh  limit,  possibly  reaching  a  limit  which  is  worse  than  that 
attained  with  the  conventional  finite  difference  approach. 

A  final  drawback  is  the  failure  of  the  variable  mesh  to 
accurately  predict  the  known  buckling  load.  If  two  good  uni¬ 
form  meshes  are  used,  the  intersection  of  the  buckling  curves 
is  the  buckling  load.  This  method  was  found  highly  useful  in 
the  last  section.  The  intersection  of  meshes  A  and  B  occurs 
at  X  =  1.9  producing  a  solution  of  P__  =  .36656EI  /L2 ,  lower 

Ci  2 

than  the  known  solution. 

For  this  beam,  changing  from  a  uniform  to  a  variable 
mesh  produced  results  inadequate  to  justify  using  a  variable 
mesh. 

The  beam  in  Figure  33  was  remodeled  to  see  if  shrinking 

the  mesh  size  around  the  discontinuity  had  any  effect.  The 

first  set  of  meshes  used  a  mesh  size  in  the  segment  enclosing 

the  discontinuity  of  h^  =  .05L.  A  second  set  of  meshes  with 

h  *  .025L  were  constructed,  with  the  same  type  of  mesh  struc- 
2 

ture  used  as  before,  i.e.  first  three  and  then  five  nodes 
around  the  discontinuity.  The  remaining  nodes  are  then 


distributed  according  to  the  length  of  the  remaining  beam  seg¬ 
ments.  The  meshes  are  similar  to  the  first  set/  but  are  more 
tightly  packed  at  the  discontinuity  and  more  loosely  packed 
away  from  the  inertia  step  change.  For  three  nodes  around  the 

step,  the  mesh  sizes  are  h  =.0875L,  h  =.025L,  and  h  =  .12917L. 

12  3 

This  is  the  A  mesh.  The  B  mesh,  with  5  nodes  about  the  dis¬ 
continuity,  has  mesh  sizes  h  =.075L,  h  =.025L,  and  h  =.1875L. 

12  3 

Checking  the  buckling  loads  at  A=1.6  (which  is  approximately  the 
optimal  wavelength  for  the  uniform  case  of  nine  nodes)  produces: 
Mesh  A  Pcr  =  .36712  Error  =  .067% 

Mesh  B  Pcr  =  .36775  Error  =  .239% 

The  buckling  load  curves  are  plotted  in  Figure  38. 

Comparing  these  results  with  the  A  and  B  meshes  from  set 
one  produces  some  interesting  results.  For  three  nodes  about 
the  discontinuity,  the  denser  mesh  of  the  second  set  is  closer 
to  the  uniform  solution  for  all  X.  However,  five  nodes  about 
the  discontinuity  favors  the  looser  mesh  of  the  first  set. 
Possibly  there  is  a  size  for  the  mesh  outside  of  the  discontin¬ 
uity  above  which  the  solution  deteriorates.  For  the  three  node 
sets,  mesh  set  two  had  higher  mesh  size  but  better  packing 
around  the  discontinuity.  Therefore  it  may  not  have  crossed 
the  limit  on  mesh  size  away  from  the  step  change  and  would 
then  produce  better  results.  On  the  other  hand,  with  the  five 
node  examples,  mesh  set  two  had  better  packing  around  the  dis¬ 
continuity  but  much  poorer  meshes  in  other  parts  of  the  beam 
than  did  the  first  case. 

Using  more  and  more  of  the  available  nodes  to  model  the 
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discontinuity  leaves  fewer  nodes  to  model  the  remaining  beam. 
This  is  the  reason  for  the  poor  performance  of  the  variable 
mesh.  The  uniform  mesh  does  a  better  job  of  modeling  the 
outer  beam  segments  while  doing  an  adequate  job  of  modeling 
the  discontinuity.  A  point  to  remember  is  that  the  disconti¬ 
nuity  needs  a  node  at  that  position.  Nodes  on  either  side 
may  create  too  high  an  inertia  gradient,  thereby  preventing 
the  buckling  curve  from  passing  through  the  critical  buckling 
load. 


Further  examples  of  the  variable  mesh  were  examined. 

The  beam  in  Figure  36  was  examined  using  three  and  five  nodes 
around  the  discontinuity,  meshes  A  and  B  respectively. 


Fig  36.  Discontinuous  Beam  P„  *  2.067EI  /L2 
3  cr  2 


Once  more,  the  uniform  nine  node  solution  was  superior 
for  all  values  of  X  except  that  small  range  where  the  variable 
mesh  curves  pass  through  the  known  solution.  The  curves  for 
Figure  36  are  shown  in  Figure  39.  The  mesh  sizes  used  were 
h  ■  .1125,  h  *  .05,  h  =  .1125  (Mesh  A)  and  h  =  .133, 


h  =  .05/  h  *  .1333  (Mesh  B) . 

2  3 

Final  checks  were  made  using  the  beams  with  exponential ly- 
and  sinusoidally-varying  inertia  from  the  previous  section. 

For  each  beam,  three  nodes  were  grouped  about  the  center  X  = 
•5L.  It  is  at  the  center  that  the  inertia  is  maximum.  Errors 
were  checked  for  two  values  of  X:  X  *  1.0  and  1.8.  The  re¬ 
sults  are  in  Table  6-1.  These  spot  checks  reveal  the  same 
pattern  as  was  previously  shown.  Namely,  that  a  variable 
mesh  produces  poorer  results  than  a  uniform  mesh  with  the 
same  number  of  nodes.  The  buckling  load  curves  for  these  two 
beams  are  available  in  Figure  40  and  Figure  41. 


TABLE  II 

Nine  Internal  Nodes 


Uniform  Mesh 
Variable  Mesh 


Sinusoid 


X  *  1.0  X  =  1.8 
.017%  -.551% 

.257%  -.663% 


Exponential 

X  =  1.0  X  =  1.8 
.449%  -.122% 

.744%  -.141% 


Summary 

For  an  equal  number  of  nodes,  a  uniform  mesh  is  generally 
much  more  accurate  than  a  variable  mesh.  The  only  place  where 
the  variable  mesh  is  more  accurate  is  the  small  X  range  where 
the  buckling  curve  passes  through  the  known  solution. 

The  variable  mesh  requires  the  buckling  wavelength 
parameter  to  be  input  before  a  buckling  load  can  be  returned. 
This  is  different  from  the  uniform  mesh  which  calculates  a 
buckling  curve  based  on  only  one  eigenvalue  matrix  solution. 


91 


ggL-tbaMJ-ac? 


Therefore,  plotting  a  variable  mesh  buckling  curve  requires 
as  many  matrix  solutions  as  points  on  the  curve,  an  extremely 
costly  method.  In  addition,  the  intersection  of  two  variable 
mesh  buckling  curves  does  not  yield  the  correct  buckling  load, 
whereas  two  uniform  mesh  curves  do.  Since  the  aim  is  to  find 
the  buckling  load  of  a  beam,  it  is  best  to  plot  two  uniform 
mesh  curves  and  find  their  intersection.  This  costs  as  much 
as  finding  two  points  on  one  variable  mesh  curve.  If  all 
that  is  desired  is  a  Pcr  for  a  choice  of  X,  the  uniform  mesh 
will  provide  a  much  better  solution  in  almost  all  cases.  The 
exception  occurs  if  X  is  chosen  near  the  optimal  X  for  the 
variable  mesh.  But  the  uniform  mesh  is  very  close  to  the 
same  solution,  so  no  real  improvement  in  accuracy  occurs. 

For  any  other  X,  a  uniform  mesh  can  provide  a  large  accuracy 
improvement.  In  fact,  for  very  large  X,  the  variable  mesh 
may  be  even  poorer  than  a  conventional  finite  difference 
solution. 
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Fig  39.  Uniform  and  variable  mesh  curves  for  a  dis¬ 
continuous  beam,  Ii/I2  =  .5,  a/L  =  . 5,  N  =  9. 
Variable  mesh  1  clusters  3  nodes  around  the 
discontinuity,  mesh  2  clusters  5.  h2  =  .05L 
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VII.  Extension  to  Composite  Materials 


Introduction 

Throughout  the  previous  sections,  attention  was  given  to 
varying  the  geometrical  properties  of  the  beam.  No  attention 
was  given  to  the  material  itself.  Material  properties  were 
assumed  isotropic,  with  the  result  that  the  beam  stiffness  E 
was  constant.  However,  the  material  stiffness  was  included 
as  a  possible  function  of  distance  along  the  beam.  Stiffness, 
in  combination  with  the  inertia,  form  a  product  known  as  the 
flexural  rigidity  El.  By  keeping  E  within  the  integration  of 
the  virtual  work  equation,  Eq  (3-30),  the  flexural  rigidity 
can  be  considered  a  function  of  X  also.  The  importance  of 
this  concept  will  be  discussed  shortly  in  con junction  with 
composite  materials. 

Composite  materials  consist  of  two  or  more  constituent 
materials  bonded  together  so  that  the  gross  material  propert¬ 
ies  are  superior  to  the  constituents.  Desirable  properties 
(high  strength,  high  stiffness,  low  weight)  are  maintained 
while  undesirable  properties  are  suppressed.  By  selecting 
the  proper  materials  and  combining  them  in  an  efficient 
geometrical  arrangement,  the  desired  gross  properties  can  be 
achieved  (35) . 

Because  composite  materials  represent  a  way  for  the 
engineer  to  design  structural  elements  to  certain  specifications 


and  still  have  control  over  other  properties  in  the  element, 
a  method  of  modeling  the  composite  beam  will  be  developed  in 
this  section. 

Theory 

A  composite  beam  is  somewhat  different  from  the  isotropic 
beams  used  in  the  previous  sections.  To  begin  with,  a  composite 
beam  can  be  formed  by  laying-up  thin  layers  of  material,  one  on 
top  of  the  other.  A  layer  of  the  beam,  called  a  lamina,  gener¬ 
ally  has  fibers  embedded  within  a  matrix.  The  fibers  provide 
most  of  the  material  strength,  with  the  matrix  serving  to  hold 
the  fibers  in  place.  As  a  lamina  is  layed-up,  the  fibers  may 
or  may  not  line  up  with  the  beam's  principal  axis.  Figure  42 
demonstrates  fibers  at  an  angle  0  to  the  beam’s  X-axis.  The 
combination  of  layers  of  material  with  different  properties 
and  orientations  from  lamina  to  lamina  requires  a  new  model 
of  the  beam. 

To  model  a  composite  beam,  two  approaches  are  available: 

1)  Laminated  beam  approach 

2)  Laminated  plate  approach 

Although  the  second  method  is  not  used  in  this  analysis,  it 
is  worth  discussing  briefly. 

Laminated  Plate  Approach 

Much  work  has  been  done  in  conjunction  with  buckling  of 
composite  plates  (28,33-36).  To  apply  the  method  to  a  beam, 
a  plate  can  be  made  very  long  and  narrow.  Two  opposite  edges 
of  the  plate  are  free  while  the  other  two  are  held  in  place. 
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Figure  42.  Typical  composite  beam  lay-up 


{ 


This  approximates  the  boundary  conditions  on  the  beam.  If 
the  free  edges  are  perpendicular  to  the  y-axis,  there  will 
be  no  moment  at  these  edges.  Because  the  plate  is  narrow. 
My  =  0  is  assumed  for  the  entire  width.  Also,  the  line  load 
along  the  plate's  x-axis  Nx  is  assumed  constant. 

The  produce  Nx  •  b,  where  b  is  the  plate  width,  approx¬ 
imates  the  axial  load  P  upon  a  beam. 

If  only  Mx  and  Nx  are  applied  to  this  narrow  plate,  a 
beam  loaded  axially  is  now  approximated.  Using  the  equations 
relating  moment,  curvature,  and  in-plane  loading  available 
in  composite  textbooks  (30),  an  expression  for  the  moment  can 
be  derived 


m  =  _L_  ifw 
x  Djj  dx2 


Sii.  N 
Dll  X 


(  7-1) 


The  coefficients  Dj ^  abd  BJ ^  are  the  bending  and  coupling 

stiffnesses  for  a  composite  laminate.  If  a  beam  of  unit  width 

is  assumed,  Nx  =  P.  For  simplicity,  the  laminated  beam  is 

assumed  symmetric  about  the  midsection.  This  simplification 

results  in  the  coupling  stiffness  Bj  being  set  equal  to  zero. 

In  addition,  the  bending  stiffness  reduces  to  D  ,  which  is 

i  i 

easily  determined  once  the  lay-up  of  the  beam  is  defined.  The 
moment  equation  has  been  reduced  to  a  result  which  is  directly 


substituted  into  the  computer  program  for  this  thesis 

d2w 

Mx  "  D i idx2 


(7-2) 


D  can  be  considered  as  an  equivalent  flexural  rigidity  El. 
A  further  derivation  is  available  in  Appendix  D. 


( 
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Results 


Very  few  results  for  the  buckling  load  of  a  composite 
column  are  available  in  literature.  Most  work  has  been  done 
in  the  area  of  plate  buckling. 

Calcote  (28)  does  provide  a  simple  theoretical  solution 
for  the  buckling  load  of  a  pinned-pinned  beam. 

d2  2 

Ncr=(Dn‘  a)?  (7“3) 

1 1 

The  buckling  stress  Ncr  is  expressed  as  load/width.  For  a 
beam  of  width  1.0,  the  total  buckling  load  p„_  =  N^„. 

Dietz  derives  a  similar  expression  for  the  buckling  load 
of  a  uniform  composite  beam  of  width  1.0  (27).  If  a  beam  has 
plies  that  are  layed-up  symmetrically  about  the  mid-section, 
it  can  be  shown  that  B  n=  0  in  Eq  (7-1).  Dietz  demonstrates 
that  for  such  a  symmetrically  laminated  beam, 

pcr  =  *2q/L’2  (7-4) 

The  (L')2  is  the  reduced  length  of  the  beam,  the  distance 

between  two  consecutive  inflection  points  of  the  deflected 

mode  shape.  L'  has  the  same  meaning  as  X  ^  did  for  a  simple 

opt 

uniform  beam.  The  optimal  wavelength  for  simple  beams  was 
found  to  correspond  to  half  the  buckled  mode's  wavelength, 
which  corresponds  to  Dietz's  L'.  Therefore,  since  X  is  known 
(Table  II  )  for  simple  cases,  Eq  (7-4)  can  be  used.  To  ana¬ 
lyze  simple  composite  beams,  the  only  change  in  solving  the 
eigenvalue  problem  is  to  replace  El  by  El  or  Dj  j  (if  symmetric). 
Therefore,  the  results  of  this  section  are  mostly  a  check  on 
calculating  material  properties. 


As  an  example,  a  composite  beam  was  made  out  of  20  plies 
of  Boron-epoxy  all  oriented  at  0°.  Material  properties  were 
obtained  from  Reference  37.  Because  the  plies  are  symmetric, 
the  buckling  load  for  a  pinned-pinned  beam  (L'  =  1)  from  Eq 
(7-4)  is 

P  *  7584. 43/L2  (lbs) 
cr 

for  a  D„  =  768.468  in- lbs.  From  a  buckling  curve,  this  load 
is  returned  for  a  X  =  1.0  (as  expected  for  a  pinned-pinned 
beam)  . 

If  the  plies  are  all  oriented  at  90°  to  the  beam's  x- 
axis,  weaker  properties  are  presented  towards  the  axial  load. 
One  would  expect  a  lower  buckling  load  because  a  boron-epoxy 
lamina  is  nine  time  less  stiff  at  90°  than  at  0°.  The  fibers 


are  not  longitudinally  aligned.  At  X  =  1.0,  a  buckling  load 
of  1028. 6/L2  (lbs)  is  returned,  matching  the  expected  solution 
from  Eq  (7-3)  or  Eq  (7-4). 

Intermixing  0°  and  90°  plies  returns  a  buckling  load  be¬ 
tween  the  two  cases  already  looked  at. 

The  buckling  load,  determined  by  knowing  X  or  find- 

opt 

ing  the  intersection  of  buckling  load  curves,  is  not  the  best 
answer.  Results  in  this  section  used  a  laminated  beam  approach 


which  assumed  =  0,  good  only  for  orientations  near  0°  and 
90°.  Results  for  other  orientations  may  not  be  as  good. 


Also,  the  effect  of  the  shear  modulus  G  has  been  ignored 


throughout  this  thesis  in  order  to  simplify  the  problem. 
Brunelle  (36)  shows  that  the  buckling  load  for  a  pinned-pinned 


uniform  beam  is 
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a  • 

4 


l 


p 

cr 


(7-5) 


where  h  is  the  beam  height.  In  this  thesis  E/G  was  assumed 
zero,  a  reasonable  assumption  for  isotropic  material.  For 
E/G  *  0,  Eq  (7-5)  returns  the  solution  discussed  in  Section 
IV.  A  composite  beam  has  a  much  higher  E/G  ratio  than  a 
solid  metal  column.  For  the  boron-epoxy  beam  with  all  layers 
at  0°,  the  ratio  E/G  is  about  40.  To  accomodate  the  trans¬ 
verse  shearing  modulus  G,  a  correction  factor  is  suggested. 
(28,33)  which  yields  the  buckling  load  for  a  pinned-pinned 


beam  as 


cr 


QTT2Er2\ 
4L2G  / 


(7-6) 


A  factor  depending  on  the  cross-section,  rj ,  is  included  as 
well  at  the  radius  of  gyration  r.  Therefore,  the  buckling 
load  returned  by  the  methods  of  this  thesis  are  going  to  be 
higher  than  they  actually  are. 

In  addition  to  the  problems  just  mentioned,  there  are 
many  others.  The  glue  holding  layers  of  the  beam  can  fail 
leading  to  gaps  in  the  beam.  Individual  fibers  within  a 
lamina  can  buckle  on  their  own  or  even  fail.  The  entire  beam 
may  fail  before  buckling  occurs.  The  solution  method  used 
in  this  thesis  does  give  a  useful  first  try  at  the  buckling 
load  of  a  composite  column  though. 


Summary 

An  approach  was  examined  with  which  a  composite  beam 
can  be  modeled.  The  laminated  beam  model  was  chosen.  This 
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model  has  the  advantage  of  direct  substitution  of  eT  for  El 
of  an  isotropic  beam.  Therefore,  isotropic  beam  results  are 
usable  for  composites.  The  solutions  obtained  matched  the 
predicted  results,  as  expected,  because  no  real  changes  are 
made  to  the  approach  of  Section  III.  Results  will  tend  to 
be  off  from  the  actual  solutions  because  of  the  laminated 
beam  approach  and  the  effect  of  the  shear  modulus  G.  How¬ 
ever,  the  solutions  achieved  can  be  useful  as  a  first  approx¬ 
imation  with  which  optimum  orientations  and  geometries  can 
be  investigated. 
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4. 


VIII.  Conclusions 

This  thesis  investigated  the  results  of  using  a  trigo¬ 
nometric  finite  difference  approach  to  solve  the  virtual  work 
equation  for  the  critical  buckling  load.  When  possible,  the 
conventional  finite  difference  solution  was  examined.  Vari¬ 
ous  values  of  X  were  used  in  the  trigonometric  approach  to 
determine  both  the  correct  solution  as  well  as  the  range  over 
which  the  trigonometric  approach  gave  more  accurate  approxi¬ 
mations  than  the  conventional  approach.  A  wide  range  of 
boundary  conditions  was  included  in  this  investigation.  In 
addition,  the  effect  of  a  variation  in  beam  inertia  and  the 

« 

use  of  a  variable  mesh  size  were  checked.  Finally,  the  appli¬ 
cation  of  the  buckling  load  technique  to  a  composite  beam 
was  made. 

Based  upon  a  recommendation  in  Hannah's  thesis,  the 
trigonometric  approach  was  applied  to  the  virtual  work  equa¬ 
tion  in  order  that  the  critical  buckling  load  could  be  cal¬ 
culated.  By  using  the  trigonometric  finite  differences,  the 
buckling  load  became  a  function  of  X  and  the  mesh  size,  as 
well  as  E,  I,  and  L.  Application  to  simple  uniform  beams 
provided  a  test  of  the  technique.  It  was  found  for  a  large 
range  of  X  that  the  trigonometric  approach  provided  more 
accuracy  than  with  conventional  finite  differences.  Increas¬ 
ing  the  number  of  nodes  also  increased  solution  accuracy. 
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The  optimum  wavelength  was  found  to  equal  half  the  buckled 
wavelength,  thereby  pinpointing  the  location  of  P  .  If  the 
buckled  shape  was  not  known  ahead  of  time,  X  =  1.4  provided 
a  fairly  accurate  solution  for  all  10  boundary  conditions. 

In  addition,  eigenvectors  produced  by  the  numerical  technique 
were  found  to  match  exactly  with  known  mode  shapes,  as  long 
as  a  node  was  positioned  at  the  point  of  maximum  displacement. 

The  effect  of  variations  in  inertia,  both  continuous  and 
discontinuous,  were  examined.  Beams  with  a  larger  middle 
cross  section,  as  in  the  case  of  the  exponential  and  sinusoi¬ 
dal  variations,  buckle  at  higher  loads  than  their  uniform 
cousins.  The  optimal  wavelength  does  not  correspond  to  any 
physical  phenomena,  such  as  the  half -buckled  wavelength.  If 
the  minimum  inertia  is  normalized  to  1.0,  the  optimal  wave¬ 
length  is  larger  than  for  the  uniform  beam.  Normalizing  the 
maximum  inertia  to  1.0  produces  an  optimal  wavelength  less 
than  the  uniform  case.  The  difference  between  the  uniform 
and  variable  inertia  beam  wavelengths  is  proportional  to  the 
size  of  the  inertia  increase.  For  the  sinusoidal  beam,  the 
inertia  did  not  change  very  rapidly.  Therefore  the  optimal 
wavelength  is  only  slightly  larger  than  for  the  uniform  beam. 
The  exponential  beam  had  larger  and  more  rapid  changes  in 
inertia.  The  choice  of  A0pt  should  be  larger  than  the  uni¬ 
form  solution  of  1.0. 

A  discontinuous  beam  was  modeled  with  maximum  inertia 
as  1.0  in  the  thesis.  Therefore,  the  optimal  wavelength 
should  be  somewhat  lower  than  the  uniform  solution,  because 
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the  variable  beam  is  less  stiff.  This  guideline  is  verified 
by  the  examples  checked. 

For  beams  with  changes  in  inertia,  the  most  effective 
buckling  load  curves  were  generated  when  a  node  was  located 
at  the  position  of  maximum  inertia  or  a  discontinuity  in  in¬ 
ertia.  By  finding  the  intersection  of  two  buckling  load 
curves  with  similar  nodal  arrangements  (i.e.,  same  points  of 
interest  modeled) ,  the  buckling  load  and  XQpt  were  establish¬ 
ed.  Failure  to  use  similar  nodal  arrangements  could  lead  to 
a  poor  solution  and  sometimes  no  solution  at  all.  The  inter¬ 
section  technique  is  better  than  guessing  at  a  wavelength 
based  upon  inertia  distribution,  if  greater  accuracy  is  re¬ 
quired. 

The  inertia  at  a  discontinuity  was  modelled  by  a  combina¬ 
tion  of  inertias  to  either  side  of  the  discontinuity.  This 
effective  inertia  maintained  continuity  in  both  slope  and 
moment.  However,  below  a  certain  value  of  effective  inertia, 
the  critical  buckling  load  can  be  reached.  The  optimal  wave¬ 
length  will  not  be  predictable  with  any  accuracy  though. 

A  variable  mesh  was  tested  in  an  effort  to  better  model 
regions  near  a  discontinuity  or  other  points  of  interest. 

While  the  buckling  curve  is  similar  in  form  to  the  uniform 
mesh  curve,  it  is  less  desirable  to  use.  The  intersection 
of  variable  mesh  curves  with  each  other  or  with  a  uniform 
mesh  curve  occurred  at  lower  buckling  loads  than  the  actual 
solution.  In  addition,  a  variable  mesh  provided  a  larger 
error  in  approximating  P  ,  compared  to  a  uniform  mesh  with 
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an  equal  number  of  nodes,  for  most  wavelengths.  Only  in  a 
small  wavelength  range  will  a  variable  mesh  be  better.  This 
could  be  useful  to  apply,  though,  if  a  uniform  mesh  cannot 
model  all  points  of  interest.  Finally,  as  more  nodes  are 
packed  near  a  discontinuity,  the  solution  accuracy  worsens. 

In  the  last  section  of  the  thesis,  a  composite  beam  was 
analyzed  using  the  laminated  narrow  plate  approach.  The 
approach  permits  the  direct  substitution  of  an  equivalent  flex¬ 
ural  rigidity  in  place  of  the  product  El.  The  thesis  then 
treats  the  composite  as  an  isotropic  beam  with  the  ability  to 
change  its  flexural  rigidity  at  each  section.  The  effects  of 
creating  a  beam  out  of  composites  does  not  directly  enter  the 
buckling  load  solution.  Instead,  the  equivalent  flexural 
rigidity  is  determined  by  using  the  composite  properties  of 
the  beam.  The  direct  substitution  mentioned  above  permits 
the  composite  beam  to  be  handled  as  easily  as  an  isotropic 
beam.  Also,  the  assumption  of  no  shear  effects  by  ignoring 
the  shear  modulus  G  produces  a  critical  buckling  load  higher 
than  in  actuality.  The  composite  beam  approach  does  provide 
a  useful  first  analysis  of  buckling  which  can  be  extended  to 
discontinuous  and  other  types  of  composite  beams  in  the  same 
manner  as  for  an  isotropic  beam. 

The  numerical  technique  applied  in  this  thesis  provided 
accurate  solutions  for  many  types  of  problems.  The  main 
point  to  be  followed  was  to  place  nodes  at  positions  of  signi¬ 
ficance  along  the  beam.  The  intersection  of  two  buckling 
load  curves  following  this  procedure  provided  the  buckling 


load  and  optimal  wavelength  for  the  truncated  Fourier  series. 
For  a  wide  range  of  wavelengths,  the  trigonometric  solutions 
are  superior  to  those  returned  by  conventional  finite  differ¬ 
ences.  If  the  intersection  of  buckling  load  curves  is  not 
used  to  choose  a  wavelength,  the  10  known  uniform  solutions 
can  act  as  guidelines  in  doing  so.  Increasing  or  decreasing 
the  wavelength  from  the  uniform  optimal  case  will  provide  a 
fairly  good  answer.  The  amount  of  increase  or  decrease  varies 
with  the  changes  in  inertia  throughout  the  beam.  Variable 
meshes  are  not  of  much  use  in  determining  the  critical  buck¬ 
ling  load,  unless  a  uniform  mesh  cannot  be  chosen  which 
models  all  the  points  of  interest.  In  this  case  a  variable 
mesh  may  be  useful  for  some  part  of  the  wavelength  range. 
Finally,  composites  can  be  modeled  by  a  simple  change  from 
El  to  the  equivalent  flexural  rigidity,  based  upon  the 
composite  properties. 


110 


Ea'SK-j*- 


Bibliography 


1.  Stein,  M.  and  Housner,  J.  Numerical  Analysis  and  Para¬ 

metric  Studies  of  the  Buckling  of  Composite  Orthotropic 
Compression  and  Shear  Panels .  NASA  Technical  Note  D-7996 
Washington:  National  Aeronautics  and  Space  Administra¬ 

tion,  October,  1975. 

2.  Stein,  M.  and  Housner,  J.  "Application  of  a  Trigonomet¬ 
ric  Finite  Difference  Procedure  to  Numerical  Analysis  of 
Compressive  and  Shear  Buckling  of  Orthotropic  Panels." 
Computers  and  Structures.  Vol.  9:17-25,  1978. 

3.  Hannah,  S.  R.  Application  of  Trigonometric  and  Conven¬ 
tional  Finite  Difference  Approximations  to  Beam  Buckling . 
Unpublished  Thesis.  Wright-Patterson  Air  Force  Base, 
Ohio:  Air  Force  Institute  of  Technology,  December  1976. 

4.  Hannah,  S.  R.  and  Palazotto,  A.  N.  "The  Incorporation 

of  Truncated  Fourier  Series  into  Finite  Difference  Approx 
imations  of  Structural  Stability  Equations."  Computer 
and  Structures.  Vol.  9:603-607,  1978. 

5.  Calcote,  L.  R.  and  Wah,  T.  Structural  Analysis  by  Finite 
Difference  Calculus.  New  York:  Van  Nostrand  Reinhold 
Co.,  1970. 

6.  Salvador!,  M.  G.  "Numerical  Computation  of  Buckling 
Loads  by  Finite  Differences."  ASCE.  Vol.  116,  Paper 
No.  2441,  1951. 

7.  Ghali,  A.  and  Neville,  A.  M.  Structural  Analysis:  A 
Unified  Classical  and  Matrix  Approach.  Scranton,  PA  : 
Intext  Educational  Publishers,  1972. 

8.  Wang,  C.  T.  Applied  Elasticity.  New  York:  McGraw-Hill 
Book  Co. ,  1953. 

9.  Hildebrand,  F.  B.  Finite-Difference  Equations  and  Simu- 
lations.  Englewood  Cliffs,  N.J.:  Prentice-Hall,  Inc. 
1968. 

10.  Wang,  Ping-Chun.  Numerical  and  Matrix  Methods  in  Struc- 
tural  Mechanics.  New  York:  John  Wiley  and  Sons,  Inc. 
1966. 

< 


111 


1 


11.  Nansen,  R.  H.  and  DiRaxnio,  H.  "Structures  for  Solar 
Power  Satellites."  Astronautics  and  Aeronautics. 

Vol.  16,  Number  10:55-59,  October  1978. 

12.  Ziegler,  H.  Principles  of  Structural  Stability. 

Waltham,  MA:  Blaisdell  Publishing  Co. ,  1968. 

13.  Michalos,  J.  and  Wilson,  E.  N.  Structural  Mechanics 
and  Analysis.  New  York:  The  MacMillian  Co.,  1965. 

14.  Girijavallabhan,  C.  V.  "Buckling  Loads  of  Non-Uniform 
Columns."  Journal  of  Structures  Division,  ASCE. 

Vol.  95:2419-2431,  1969. 

15.  Abbassi,  M.  M.  "Buckling  of  Struts  of  Variable  Bending 
Rigidity."  Journal  of  Applied  Mechanics,  Transactions 
of  ASME.  Vol.  25:537-540,  December  1958. 

16.  Abbassi,  M.  M.  "The  Second  Approximation  for  Buckling 
Loads  of  Tapered  Struts."  Journal  of  Applied  Mechanics, 
Transactions  of  ASME.  Vol.  28:211-212,  March  1960. 

17.  Sinha,  S.  C.  and  Chow,  C.  C.  "Approximate  Eigenvalues 
for  Systems  with  Variable  Parameters."  Journal  of 
Applied  Mechanics.  Vol.  46:203-205,  March  1979. 

18.  Almroth,  B.  O.  and  Brush,  D.  0.  Buckling  of  Bars, 

Plates,  Shells.  New  York:  McGraw-Hill  Book  Co.,  1965. 

19.  Popov,  E.  P.  Introduction  to  Mechanics  of  Solids. 

Englewood  Cliffs,  N.J. :  Prentice-Hall,  Inc.,  1968. 

20.  Timoshenko,  S.  P.,  and  Young,  D.  H.  Theory  of  Structures . 

New  York:  McGraw-Hill  Book  Co.,  1965. 

21.  Willems,  N.  and  Lucas,  W.  M.  Matrix  Analysis  for  Struc¬ 
tural  Engineers.  Englewood  Cliffs,  N.J.:  Prentice-Hall , 
Inc.,  1968. 

22.  Timoshenko,  W.  P.  and  Gere,  J.  M.  Theory  of  Elastic 
Stability.  New  York:  McGraw-Hill  Book  Co.,  1961'. 

23.  Shaker,  J.  Effect  of  Axial  Load  on  Mode  Shapes  and 
Frequencies  of  Beams.  NASA  Technical  Note  D-8109, 

Washington:  National  Aeronautics  and  Space  Adminis¬ 
tration,  November  1975. 

24.  O’Rourke,  M.  and  Zebrowski,  T.  "Buckling  Load  for  Non- 
Uniform  Columns."  Computers  and  Structures.  Vol.  7:717-720, 
1977. 

25.  Acton,  F.  S.  Numerical  Methods  that  Work.  New  York: 

Harper  and  Row  Publishers,  Inc.,  19 V0. 


112 


26.  Franklin,  J.  N.  Matrix  Theory .  Englewood  Cliffs,  N.J.: 
Prentice-Hall  Inc.,  1968. 

27.  Dietz,  A.  G.  H.  (ed)  Composite  Engineering  Laminates. 
Cambridge,  MA:  MIT  Press,  1969. 

28.  Calcote,  L.  R.  The  Analysis  of  Laminated  Composite 
Structures.  New  York:  Van  Nostrand  Reinhold  Co.,  1969. 

29.  Hopkins,  R.  B.  Design  Analysis  of  Shafts  and  Beams. 

New  York:  McGraw-Hill  Book  Co.,- T9 70. 

30.  Jones,  R.  M.  Mechanics  of  Composite  Materials. 

Washington,  D.C. :  Scripts  Book  Co.,  1975. 

31.  Moler,  C.  B. ,  and  Stewart,  C.  W.  "An  Algorithm  for 
Generalized  Matrix  Eigenvalue  Problems."  SIAM 
Journal  of  Numerical  Analysis,  10,  1973. 

32.  Almroth,  B.  O.,  and  Brogan,  F.  A.  "Numerical  Analysis 
of  Structures  (Discretization  Procedures) . "  Lockheed 
Missiles  and  Space  Company,  Inc. :  Technical  Report  # 
LMSC-D556462 ,  January  1977. 

33.  Broutman,  L.  J.,  and  Krock,  R.  H.  Composite  Materials. 
Vol.  7  Structural  Design  and  Analysis  -  Part  I.  Ed.  by 
C.  C.  Chamis.  New  York:  Academic  Press,  Inc.,  1975. 

34.  Tsai,  S.  W.  Analyses  of  Composite  Structures.  NASA  Con¬ 
tract  Report  CR-620,  Washington:  National  Aeronautics 
and  Space  Administration,  November,  1966. 

35.  Tsai,  S.  W.  Structural  Behavior  of  Composite  Materials. 
NASA  Contract  Report  CR-71,  Washington:  National  Aero- 
nautics  and  Space  Administration,  July,  1964. 

36.  Brunell,  E.  J.  "The  Statics  and  Dynamics  of  a  Transversely 
Isotropic  Timoshenko  Beam."  Journal  of  Composite  Materi¬ 
als.  Vol.  4:405-407,  1970. 

37.  Mandell,  J.  R.  An  Experimental  Investigation  of  the 
Buckling  of  Anisotropic  Fiber  Reinforced  Plastic  Plates . 
Air  Force  Materials  Laboratory;  AFML-TR-68-281 ,  October 
1968;  Wright-Patterson  Air  Force  Base,  Ohio. 


113 


0i'i* 


Also  from  beam  theory,  the  stress  and  moment  are  related  by 

a  =  My/I  and  the  inertia  I  =  /  y2dy  dz.  Using  these  defini- 

A 

tions,  the  strain  energy  of  the  beam  reduces  to 


U i 

U-  J  2EI 


dx 


(A-3) 


The  moment  in  the  beam  can  be  expressed  in  terms  of  the  de¬ 
flection  as  M  a  EIw" .  The  strain  energy  reduces  further  to 
the  familiar  form 

U  =  dx  (A-4) 

The  internal  virtual  work  is  the  variation  of  strain  energy, 
expressed  as 


<5U  = 


d2w  d2 6w 
dx2  dx2 


dx 


(A-5) 


In  addition  to  the  internal  work,  external  work  is  produced 
by  the  axial  load  as  it  moves  through  the  displacement  du. 


The  external  work  Vf  is 

e 


dW  =  -JjP  du 
e 


(A- 6) 


( 


The  displacement  du  needs  to  be  expressed  in  terms  of  the 
vertical  displacement  w  before  the  external  work  expression 
is  useful.  Figure  A2  illustrates  a  segment  of  the  deflected 
beam  segment  AB'  compared  with  the  undisplaced  segment  AB. 

An  expression  for  AB'  is 

1* 


AB 


'  «  ( 


<a*>‘- (gdxj 


(A-7) 


This  expression  can  be  approximated,  through  the  use  of  the 
binomial  expansion,  as 


AB' 


dx 


/jwV 

1  _  Mr*)  +  ••• 

\dx  J 


(A-8) 
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Fig  A2 .  Deflected  Beam  Segment  Shape 


The  displacement  du  equals  the  segment  B'B,  which  in  turn 

equal,  s  /dwy 

AB  -  AB'  =  -  )  dx 

The  total  horizontal  displacement  of  the  beam  is  found  by 

integrating  Equation  9  along  the  length  of  the  beam 


(A-9) 


u  = 


(A-10) 


The  total  external  work  done  by  the  axial  load  is  found  by 
integrating  Eq  (A6)  and  substituting  Eq  (AlO) 

dw\2 


we  » 


/( 


dx/ 


dx 


(A-ll) 


The  virtual  work  of  the  external  force  is 


6W  =  P 
e 


/dw  d6w 
dx  dx 


dx  (A-12) 

For  equilibrium,  the  variation  in  total  potential 


energy  v  must  be  zero,  where  ir  =  U  -  Wfi.  This  gives  the 


virtual  work  equation 
and 


6  U  =  6We 


(A-13 


L 

/  El 
0  dx2  dx2 


d2w  d26w 


dx  =  P  /  dx 

0  dx  dx 


(A-14 


APPENDIX  B 


Trigonometric  Finite  Difference  Expression  For  the  Second 
Derivative  of  Deflection  with  Variable  Mesh 

When  a  uniform  mesh  was  used,  the  expression  for  the 
second  derivative  at  node  i  was  developed  in  Section  III  as 

Wi  =  |z(Wi+l  "  2Wi  +  wi-l}  (B-l) 

This  derivative  applies  only  when  the  mesh  size  is  constant 

across  the  node-  In  Section  VI,  a  variable  mesh  is  used 
where  the  second  derivative  is  required  at  a  discontinuity 
in  mesh  size.  A  new  expression  is  required  for  the  second 
derivative  at  the  change  in  mesh  size.  Figure  B1  shows  the 
mesh  at  the  discontinuity,  with  mesh  sizes  h'  and  h". 


Fig  Bl.  Variable  Mesh  Size  Arrangement 


The  trigonometric  expression  used  to  generate  finite 
differences  is 

,  tMX-X.)  / 1  \  2  r  "(X-X.n 


Applying  the  above  expression  at  X  =  X,  +  h"  and  X  =  X^  -  h' 
yields  the  equations 

;’)WV  (B— 3 ) 


wi+l 

"  Wi 

X  .  Trh'*,, 

+  *sln“  Wi 

♦  SO 

irh' 

L  -  cos— 

w.  , 

=  W. 

_  Asin*f  W! 

ir  X  i 

irh 

L  -  cos— 

i-1 

l 

*  V 

X 

The  needed  expression  for  W7  can  be  obtained  with  the  follow- 

Trh" 

mg  procedure.*  Multiply  by  sin—  ,  multiply  W.  by 

.  irh 1 

sm— ^  and  add  the  two  expressions.  Using  the  expression 


.  irh  ~  .  »irh 

1  -  cos—  =  2sin 


reduces  the  second  derivative  to  the  desired  result, 
.irh 


sin- 


W7  = 

l 


L’wi+1  -  (si 


sin—  +  sin 


in—"  +  sin— "w^i 


,/xV  .  TTh'  .  2irh" 

[sinT  Sin  21 


.  Trh"  .  2  irh 
+  sin—  sin22X 


I  * 


(B-5) 


Equation  (B5)  is  the  expression  required 
ative  at  a  node  where  the  mesh  changes. 
(B5)  reduces  to  Eq  (Bl)  where 


Z  21  .  Trh 
h  TT  Sln2TT 


for  the  second  deriv- 
If  h'  =  h"  =  h,  Eq 
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APPENDIX  C 


Effective  Inertia  at  a  Sudden  Change  in  Cross-Section 


For  a  beam  with  a  discontinuity  at  node  i,  as  shown 
in  Figure  Cl, ‘what  is  the  value  of  the  inertia  at  that  node? 


Fig  Cl.  Step  Change  in  Inertia 
The  inertia  of  the  section  left  of  the  discontinuity  is  IL 
and  to  the  right,  IR.  With  these  values  and  the  deflection 
curve  in  Figure  C2,  the  effective  inertia  can  be  calculated. 
The  curves  in  Figure  C2  are  beam  displacements  from  node 
i  -  1  to  node  i  +  1.  Curve  ABC  is  the  true  deflection,  with 
displacements  W W^,  W^+1,  respectively.  If  curve  A-B  is 
extended  to  point  C',  a  fictitious  displacement  is  gen¬ 

erated  for  node  i  +  1.  Likewise,  extending  curve  C-B  to  point 
A'  produces  a  fictitious  displacement  ^  at  node  i  -  1. 
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Fig  C2 .  Fictitious  and  Actual  Node  Displacements 
For  compatibility,  the  slope  of  the  two  curves  at  i 
must  be  equal.  Using  a  full-station  central  difference  for 
slope 

Wi  -  2H<Ki+l  -  wi-l>  (C-l) 

yields  the  continuity  equation 

Wi+1  '  Wi-1  =  wi+l  "  wi-l  (C-2) 


Also,  for  equilibrium,  the  bending  moments  either  side  of 
node  i  must  be  equal.  The  moment  is  described  by  M  =  El^pr 
Using  a  full-station  central  difference  to  equate  moments 
produces 


eil  f 

Mf  =  -^(W?+1-2Wi  +  W^) 
EX 

mi  -  ^<wi+i-3wi +  "Li* 


(C-3) 
(C-4 ) 


Expressing  in  terms  of 


from  Eq  (C2)  and  inserting 


into  Eq  (C3)  and  Eq  (C4)  produces  two  equations  for  in 
f 

terms  of  w^+^*  Eliminating  the  fictitious  displacement  from 
both  equations  produces 


M1  h2(Ir  +  In) 


L  '  J'R 

The  general  form  for  the  moment  at  node  i  is 


Ki+i  -  2wi  +  Wi-J  (c-5) 


EIeff  r  1 

Mi  -  17-  LWi+l  "  2Wi  +  Wi-lJ 


(C-6) 


By  comparing  the  two  moment  equations,  the  effective  inertia 
at  a  discontinuity  is 


I 


eff 


2ilxr 

XL  +  XR 


(C-7) 
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APPENDIX  D 


Flexural  Rigidity  of  a  Composite  Beam 


The  laminated  beam  approach  is  used  to  model  the  flex¬ 
ural  rigidity  El  of  a  composite  beam. 


From  this  approach 


and 


1  1  i  „i 

az  °xy  ~  ayz  ~  My 


a1  =  qH1 

X  i!  X 


(D-l) 


(D-2 ) 


The  superscript  i  is  the  layer  number  for  a  ply  in  the  com¬ 
posite.  The  assumptions  in  Eg  (Dl)  yield  the  stress  equation 
(D2).  Stress  in  each  layer  is  modeled  by  Hookes'  Law.  The 
strain  can  be  expressed  in  terms  of  distance  from  the  mid¬ 
section  and  the  curvature  as  Z(d2w/dx2).  Substituting  into 
Eg  (D2)  gives 


=QmZ 


(D-3) 


The  moment  in  a  simple  isotropic  beam  is  given  by  the  formula 


M  =  El 


d2w 

dx2 


(D-4 ) 


For  a  composite  beam,  the  moment  is  found  by  integrat¬ 
ing  the  stress  and  its  moment  arm  throughout  the  thickness. 

t/2  ± 

M  =  b  /  <j^Z  dz  (D-5) 

-t/2  X 

The  composite  beam  is  of  constant  width  b  and  thickness  t. 
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1 


Substituting  for  a*  from  Eq  (D3)  produces 


M  = 


t/2_  i 
/  Q  Z2dz 
-t/2  11 


\  d2w 

ydx2 


(D-6) 


Comparing  Eq  (D4)  and  Eq  (D6)  shows  that  the  same  moment- 
curvature  behavior  results  for  a  composite  beam  if  an  equi¬ 
valent  flexural  rigidity 


_  t/2_  . 

El  =  b  /  Q  xZ2dz  (D-7) 

-t/2  11 


is  substituted  into  Eq  (D6).  The  moment  equation  becomes 


M 


(D-8) 


Because  the  eigenvalue  equation  of  this  thesis  was  derived 
using  Eq  (D4),  a  simple  substitution  of  El  for  El  permits 
modeling  of  composite  beams. 

Using  this  equation,  Dietz  (28)  reduces  the  equivalent 
flexural  rigidity  to  the  form 


The  values  of  D  ,  B  ,  A  are  dependent  upon  the  thickness, 
111111 

location,  orientation,  and  material  properties  of  each  layer. 
The  equations  describing  1}  ,  ^  ,  and  ^  t  are  available  in 
composite  textbooks  (30) . 


( 
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